Atomistic modelling of the immiscible Fe–Bi system from a constructed bond order potential

Z B Liang; Y C Jiang; X Gong; H R Gong

doi:10.1088/1361-648X/ac2e8e

1. Introduction

Iron–bismuth (Fe–Bi) is an immiscible binary system with a positive heat of formation (+18 kJ mol⁻¹) [1] and has raised a lot of research interests during the past decades. On the one hand, Fe–Bi nanoparticles and multilayers possess unique magnetic properties and have been potentially used in some industrial and biomedical fields such as computer tomography and magnetic resonance imaging [2–7]. Interestingly, metastable rhombohedral and body-centered-cubic (BCC) solid solutions as well as amorphous structures of Fe–Bi have been experimentally fabricated by ion beam mixing [8], magnetron sputtering [3, 4, 7], electron-beam evaporation [2, 9, 10], mechanical alloying [11], etc. On the other hand, the interaction between Fe and Bi seems important in nuclear power systems, as Bi and its alloys are a first-choice candidate of spallation target and coolant, which induce the corrosion and liquid metal embrittlement (LME) of stainless steels and affect the safety operation of the nuclear reactors [12–14].

It is well known that molecular dynamics (MD) simulation based on a realistic n-body potential is a capable way to predict various structures and properties of a binary system. In this respect, two Fe–Bi potentials have been already constructed within the framework of Lennard-Jones (LJ) and have been used to study the corrosion behavior of the Fe–Bi system [15, 16]. Moreover, a Bi–Bi potential has been derived under the formalism of embedded atom model (EAM) to simulate various properties of liquid Bi [17, 18]. However, we will show later that the above three potentials are failed to describe the ground-state rhombohedral structure (A7, space group No. 166, R $\bar{3}$ m) of Bi [15–18] and are unable to predict the structural stability of metastable phases formed in the immiscible Fe–Bi system [15, 16].

Bond-order-potential (BOP) has been well regarded as a high-fidelity n-body potential based on molecular orbital theory and quantum mechanical tight binding models [19, 20]. As an angular dependent potential, BOP is especially appropriate for strong covalent-bonding materials such as Si, C, GaAs [21–24] and is also suitable for transition metals and alloys [19, 20, 25]. To our knowledge, however, there is not any BOP potential of Bi published so far in the literature. Therefore, a new n-body Fe–Bi potential is first constructed in the present study under the formalism of BOP [19, 20, 26–29]. Based on this Fe–Bi potential, MD simulations are performed to reveal the ground-state properties of the A7 structure of Bi, the structural stability of metastable phases of Fe–Bi, and the effects of Bi on mechanical properties of BCC Fe. It will be shown that the newly constructed Fe–Bi potential could correctly predict structural parameters and various properties of the ground-state A7 structure of Bi. The derived results from the present Fe–Bi potential simulations will be also discussed and compared with similar observations from the other potentials, experiments, and thermodynamical models in the literature [1, 15–18].

2. Construction of Fe–Bi bond-order potential

As related before, two LJ potentials and an EAM potential have been constructed for Bi and Fe–Bi [15–18], while it will be shown later that they are unstable to reproduce the ground-state rhombohedral structure (A7, space group No. 166, R $\bar{3}$ m) of Bi as well as structural stability of Fe–Bi. In the present study, the Bi–Bi, Fe–Fe, and Fe–Bi potentials are therefore constructed based on the framework of BOP by Pettifor et al [26–29], and are compared extensively with the Fe–Bi potentials already in the literature.

According to the BOP formalism, the bond energy is the product of the bond-order and bond integral, and the total energy of a whole system with N atoms is expressed as [20]:

$\begin{align}E& =\frac{1}{2}\sum\limits _{i=1}^{N}\sum\limits _{j={i}_{1}}^{{i}_{N}}{\varnothing }_{ij}\left({r}_{ij}\right)-\sum\limits _{i=1}^{N}\sum\limits _{j={i}_{1}}^{{i}_{N}}{\beta }_{\sigma ,ij}\left({r}_{ij}\right)\cdot {{\Theta}}_{\sigma ,ij}\\ & \quad -\sum\limits _{i=1}^{N}\sum\limits _{j={i}_{1}}^{{i}_{N}}{\beta }_{\pi ,ij}\left({r}_{ij}\right)\cdot {{\Theta}}_{\pi ,ij},\end{align} \tag{ 1 }$

where r_ij is the distance between atoms i and j, ∅_ij, β_σ,ij, and β_π,ij are three pairwise functions, Θ_σ,ij and Θ_π,ij are two variables which describe local circumstances around atoms i and j, and the list j = i₁, i₂, ..., i_N represents neighbors of atom i. The first term of the right-hand side of equation (1) represent the overlap repulsion between a pair of ion cores [30]. β_σ,ij and β_π,ij are σ and π bond integrals, respectively, while Θ_σ,ij and Θ_π,ij are σ and π bond orders, respectively. The three pairwise funcitons have similar forms as:

$\begin{equation}{\varnothing }_{ij}\left({r}_{ij}\right)={\varnothing }_{0,ij}\cdot {f}_{ ij}{\left({r}_{ij}\right)}^{{m}_{ij}}\cdot {f}_{ \mathrm{c},ij}\left({r}_{ij}\right),\end{equation} \tag{ 2 }$

$\begin{equation}{\beta }_{\sigma ,ij}\left({r}_{ij}\right)={\beta }_{\sigma ,0,ij}\cdot {f}_{ij}{\left({r}_{ij}\right)}^{{n}_{ij}}\cdot {f}_{ \mathrm{c},ij}\left({r}_{ij}\right),\end{equation} \tag{ 3 }$

$\begin{equation}{\beta }_{\pi ,ij}\left({r}_{ij}\right)={\beta }_{\pi ,0,ij}\cdot {f}_{ij}{\left({r}_{ij}\right)}^{{n}_{ij}}\cdot {f}_{\mathrm{c},ij}\left({r}_{ij}\right),\end{equation} \tag{ 4 }$

where ∅_0,ij, β_σ,0,ij, β_π,0,ij, m_ij, and n_ij are five pair-dependent parameters needed to be fitted. ∅_0,ij, β_σ,0,ij, β_π,0,ij are three prefactors of repulsive energy, σ bond integral, and π bond integral, receptively. ${f}_{ ij}\left({r}_{ij}\right)$ is a Goodwin–Skinner–Pettifor (GSP) radial function written as [21]:

$\begin{equation}{f}_{ij}\left({r}_{ij}\right)=\frac{{r}_{0,ij}}{{r}_{ij}}\enspace \mathrm{exp}\left[{\left(\frac{{r}_{0,ij}}{{r}_{\mathrm{c},ij}}\right)}^{{n}_{\mathrm{c},ij}}-{\left(\frac{{r}_{ij}}{{r}_{\mathrm{c},ij}}\right)}^{{n}_{\mathrm{c},ij}}\right],\end{equation} \tag{ 5 }$

in which r_0,ij, r_c,ij, and n_c,ij are three pair-dependent parameters, and r_0,ij and r_c,ij are GSP reference radius and GSP characteristic radius, respectively. ${f}_{\mathrm{c},ij}\left({r}_{ij}\right)$ is a cutoff function defined as:

$\begin{equation}{f}_{\mathrm{c},ij}\left({r}_{ij}\right)=\begin{cases}\frac{\mathrm{exp}\left(-{\alpha }_{ij}\cdot {{r}_{ij}}^{{\gamma }_{ij}}\right)-\mathrm{exp}\left(-{\alpha }_{ij}\cdot {{r}_{\mathrm{c}\mathrm{u}\mathrm{t},ij}}^{{\gamma }_{ij}}\right)}{\mathrm{exp}\left(-{\alpha }_{ij}\cdot {{r}_{1,ij}}^{{\gamma }_{ij}}\right)-\mathrm{exp}\left(-{\alpha }_{ij}\cdot {{r}_{\mathrm{c}\mathrm{u}\mathrm{t},ij}}^{{\gamma }_{ij}}\right)},\quad {r}_{ij}< {r}_{\mathrm{c}\mathrm{u}\mathrm{t},ij}\quad \\ 0,\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad {r}_{ij}\geqslant {r}_{\mathrm{c}\mathrm{u}\mathrm{t},ij}\quad \end{cases},\end{equation} \tag{ 6 }$

where r_1,ij and r_cut,ij are two pair-dependent parameters to represent cutoff start radius and cutoff radius, respectively, and γ_ij and α_ij can be obtained from ${\gamma }_{ij}=\frac{\mathrm{ln}\left[\mathrm{ln}\left(0.99\right)/\mathrm{ln}\left(0.01\right)\right]}{\mathrm{ln}\left(\frac{{r}_{1,ij}}{{r}_{\mathrm{c}\mathrm{u}\mathrm{t},ij}}\right)}$ and ${\alpha }_{ij}=-\frac{\mathrm{ln}\left(0.99\right)}{{\left({r}_{1,ij}\right)}^{{\gamma }_{ij}}}$ .

The variable, Θ_σ,ij, is defined as

$\begin{align}{{\Theta}}_{\sigma ,ij}& ={{\Theta}}_{s,ij}({{\Theta}}_{\sigma ,ij}^{(1 /2)},{f}_{\sigma ,ij})\cdot \\ & \quad \times [1-({f}_{\sigma ,ij}-\frac{1}{2})\cdot {k}_{\sigma ,ij}\\ & \quad \times \frac{{\beta }_{\sigma ,ij}^{2}({r}_{ij})\cdot {\mathfrak{R}}_{3\sigma ,ij}}{{\beta }_{\sigma ,ij}^{2}({r}_{ij})+\frac{{\beta }_{\sigma ,ij}^{2}({r}_{ij})\cdot {{\Phi}}_{2\sigma }^{i}+{\beta }_{\sigma ,ij}^{2}({r}_{ij})\cdot {{\Phi}}_{2\sigma }^{j}}{2}+{\zeta }_{2}}],\end{align} \tag{ 7 }$

where ${{\Theta}}_{\sigma ,ij}^{\left(1 /2\right)}$ is the two-level σ bond order for the BOP approximation with a half-full valence shell, ${{\Phi}}_{2\sigma }^{i}$ and ${{\Phi}}_{2\sigma }^{j}$ are local variables arising from the 2-hops, and ${\mathfrak{R}}_{3\sigma ,ij}$ is a local variable arising from the three-member ring hops [31]. ζ₂ and ζ₁, ζ₃, ζ₄ in the following equations (10) and (15) are global parameters, f_σ,ij and k_σ,ij are two pair-dependent parameters to be fitted. ${{\Theta}}_{s,ij}\left({{\Theta}}_{\sigma ,ij}^{\left(1 /2\right)},{f}_{\sigma ,ij}\right)$ is a function with continuous derivatives, and is expressed as the Ward's solution [31]:

$\begin{align}& {{\Theta}}_{s,ij}\left({{\Theta}}_{\sigma ,ij}^{\left(1 /2\right)},{f}_{\sigma ,ij}\right)\\ & \qquad =\frac{{{\Theta}}_{0}+{{\Theta}}_{1}+S\cdot {{\Theta}}_{\sigma ,ij}^{\left(1 /2\right)}-\sqrt{{\left({{\Theta}}_{0}+{{\Theta}}_{1}+S\cdot {{\Theta}}_{\sigma ,ij}^{\left(1 /2\right)}\right)}^{2}-4\left(-\varepsilon \sqrt{1+{S}^{2}}+{{\Theta}}_{0}+{{\Theta}}_{1}+S\cdot {{\Theta}}_{\sigma ,ij}^{\left(1 /2\right)}\right)}}{2},\end{align} \tag{ 8 }$

where

$\begin{equation*} \varepsilon =1{0}^{-10}\end{equation*}$

$\begin{align*}& {{\Theta}}_{0}=15.737\enspace 980\cdot \\ & {\left(\frac{1}{2}-\left\vert {f}_{\sigma ,ij}-\frac{1}{2}\right\vert \right)}^{1.137\enspace 622}\cdot {\left\vert {f}_{\sigma ,ij}-\frac{1}{2}\right\vert }^{2.087\enspace 779}\end{align*}$

$\begin{align*}S& =1.033\enspace 201\cdot \left\{1-\mathrm{exp}\left[-22.180\enspace 680\right.\right.\\ & \quad \left.\left.\cdot {\left(\frac{1}{2}-\left\vert {f}_{\sigma ,ij}-\frac{1}{2}\right\vert \right)}^{2.689\enspace 731}\right]\right\} \end{align*}$

$\begin{equation} {{\Theta}}_{1}=2\cdot \left(\frac{1}{2}-\left\vert {f}_{\sigma ,ij}-\frac{1}{2}\right\vert \right).\end{equation} \tag{ 9 }$

The variable ${{\Theta}}_{\sigma ,ij}^{\left(1 /2\right)}$ in equations (7) and (8) is expressed as:

$\begin{equation}{{\Theta}}_{\sigma ,ij}^{\left(1 /2\right)}=\frac{{\beta }_{\sigma ,ij}\left({r}_{ij}\right)}{\sqrt{{\beta }_{\sigma ,ij}^{2}\left({r}_{ij}\right)+{c}_{\sigma ,ij}\cdot \left[{\beta }_{\sigma ,ij}^{2}\left({r}_{ij}\right)\cdot {{\Phi}}_{2\sigma }^{i}+{\beta }_{\sigma ,ij}^{2}\left({r}_{ij}\right)\cdot {{\Phi}}_{2\sigma }^{j}\right]+{\zeta }_{1}}},\end{equation} \tag{ 10 }$

where c_σ,ij is a pair-dependent parameter. The variable ${{\Phi}}_{2\sigma }^{i}$ in equations (7) and (10) has the following form:

$\begin{equation}{\beta }_{\sigma ,ij}^{2}\left({r}_{ij}\right)\cdot {{\Phi}}_{2\sigma }^{i}=\sum\limits _{\begin{matrix}\hfill k={i}_{1}\hfill \\ \hfill k\ne j\hfill \end{matrix}}^{{i}_{N}}{g}_{\sigma ,jik}^{2}\left({\theta }_{jik}\right)\cdot {\beta }_{\sigma ,ik}^{2}\left({r}_{ik}\right),\end{equation} \tag{ 11 }$

where θ_jik is the bond angle of atoms j–i–k, and g_σ,jik(θ_jik) is the angular function using the Ward's quadratic spline formula [20] as

$\begin{equation}{g}_{\sigma ,jik}\left({\theta }_{jik}\right)={g}_{0}+{g}_{1}\enspace \mathrm{cos}\left({\theta }_{jik}\right)+{g}_{2}\enspace \mathrm{c}\mathrm{o}{\mathrm{s}}^{2}\left({\theta }_{jik}\right),\end{equation} \tag{ 12 }$

where

$\begin{equation*}{g}_{0}=\frac{\left({b}_{\sigma ,jik}-{g}_{0,jik}\right)\cdot {u}_{\sigma ,jik}^{2}-\left({g}_{0,jik}+{b}_{\sigma ,jik}\right)\cdot {u}_{\sigma ,jik}}{2\cdot \left(1-{u}_{\sigma ,jik}^{2}\right)}\end{equation*}$

$\begin{equation*} {g}_{1}=\frac{{g}_{0,jik}+{b}_{\sigma ,jik}}{2}\end{equation*}$

$\begin{equation} {g}_{2}=\frac{{g}_{0,jik}-{b}_{\sigma ,jik}+\left({g}_{0,jik}+{b}_{\sigma ,jik}\right)\cdot {u}_{\sigma ,jik}}{2\cdot \left(1-{u}_{\sigma ,jik}^{2}\right)},\end{equation} \tag{ 13 }$

in which g_0,jik, b_σ,jik and u_σ,jik are three three-body-dependent (jik) parameters to be fitted. It should be noted that the variable ${{\Phi}}_{2\sigma }^{j}$ in equations (7) and (10) is similar to ${{\Phi}}_{2\sigma }^{i}$ with the bond angle of atoms i–j–k.

The term ${\beta }_{\sigma ,ij}^{2}\left({r}_{ij}\right)\cdot {\mathfrak{R}}_{3\sigma ,ij}$ in equation (7) is calculated as:

$\begin{align}{\beta }_{\sigma ,ij}^{2}\left({r}_{ij}\right)\cdot {\mathfrak{R}}_{3\sigma ,ij}& =\sum\limits _{\begin{matrix}\hfill k={i}_{1}\hfill \\ \hfill k,j=n\hfill \end{matrix}}^{{i}_{N}}{g}_{\sigma }\left({\theta }_{jik}\right)\cdot {g}_{\sigma }\left({\theta }_{ijk}\right)\cdot {g}_{\sigma }\left({\theta }_{ikj}\right)\\ & \quad \cdot {\beta }_{\sigma ,ik}\left({r}_{ik}\right)\cdot {\beta }_{\sigma ,jk}\left({r}_{jk}\right),\end{align} \tag{ 14 }$

where k, j = n means that atoms k and j are neighbors.

The variable Θ_π,ij in equation (1) is expressed as

$\begin{align}{{\Theta}}_{\pi ,ij}& =\frac{{a}_{\pi ,ij}\cdot {\beta }_{\pi ,ij}\left({r}_{ij}\right)}{\sqrt{{\beta }_{\pi ,ij}^{2}\left(rij\right)+{c}_{\pi ,ij}\cdot \left(\frac{{\beta }_{\pi ,ij}^{2}\left({r}_{ij}\right)\cdot {{\Phi}}_{2\pi ,ij}^{i}+{\beta }_{\pi ,ij}^{2}\left({r}_{ij}\right)\cdot {{\Phi}}_{2\pi ,ij}^{j}}{2}-\sqrt{{\beta }_{\pi ,ij}^{4}\left({r}_{ij}\right)\cdot {{\Phi}}_{4\pi ,ij}+{\zeta }_{3}}+\sqrt{{\zeta }_{3}}\right)+{\zeta }_{4}}}\\ & \quad +\frac{{a}_{\pi ,ij}\cdot {\beta }_{\pi ,ij}\left({r}_{ij}\right)}{\sqrt{{\beta }_{\pi ,ij}^{2}\left(rij\right)+{c}_{\pi ,ij}\cdot \left(\frac{{\beta }_{\pi ,ij}^{2}\left({r}_{ij}\right)\cdot {{\Phi}}_{2\pi ,ij}^{i}+{\beta }_{\pi ,ij}^{2}\left({r}_{ij}\right)\cdot {{\Phi}}_{2\pi ,ij}^{j}}{2}+\sqrt{{\beta }_{\pi ,ij}^{4}\left({r}_{ij}\right)\cdot {{\Phi}}_{4\pi ,ij}+{\zeta }_{3}}\right)+{\zeta }_{4}}},\end{align} \tag{ 15 }$

where a_π,ij and c_π,ij are two pair-dependent parameters.

The term ${\beta }_{\pi ,ij}^{2}\left({r}_{ij}\right)\cdot {{\Phi}}_{2\pi ,ij}^{i}$ is calculated as:

$\begin{align}{\beta }_{\pi ,ij}^{2}\left({r}_{ij}\right)\cdot {{\Phi}}_{2\pi ,ij}^{i}& =\sum\limits _{\begin{matrix}\hfill k={i}_{1}\hfill \\ \hfill k\ne j\hfill \end{matrix}}^{{i}_{N}}\left[{p}_{\pi ,i}\cdot {\beta }_{\sigma ,ik}^{2}\left({r}_{ik}\right)\cdot \mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\left({\theta }_{jik}\right)\right.\\ & \quad \left.+\left(1+\mathrm{c}\mathrm{o}{\mathrm{s}}^{2}\left({\theta }_{jik}\right)\right)\cdot {\beta }_{\pi ,ik}^{2}\left({r}_{ik}\right)\right],\end{align} \tag{ 16 }$

where p_π,i is an element-dependent parameter.

The term ${\beta }_{\pi ,ij}^{4}\left({r}_{ij}\right)\cdot {{\Phi}}_{4\pi ,ij}$ in equation (15) can be written as:

$\begin{equation}\begin{aligned}{\beta }_{\pi ,ij}^{4}\left({r}_{ij}\right)\cdot {{\Phi}}_{4\pi ,ij}& =\frac{1}{4}\sum\limits _{\begin{matrix}\hfill k={i}_{1}\hfill \\ \hfill k\ne j\hfill \end{matrix}}^{{i}_{N}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{4}\left({\theta }_{jik}\right)\cdot {\hat{\beta }}_{ik}^{4}\left({r}_{ik}\right)\\ & \quad +\frac{1}{4}\sum\limits _{\begin{matrix}\hfill k={j}_{1}\hfill \\ \hfill k\ne i\hfill \end{matrix}}^{{j}_{N}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{4}\left({\theta }_{ijk}\right)\cdot {\hat{\beta }}_{jk}^{4}\left({r}_{jk}\right)\\ & \quad +\frac{1}{2}\sum\limits _{\begin{matrix}\hfill k={i}_{1}\hfill \\ \hfill k\ne j\hfill \end{matrix}}^{{i}_{N}}\sum\limits _{\begin{matrix}\hfill {k}^{\prime }=k + 1\hfill \\ \hfill { k}^{\prime }\ne j\hfill \end{matrix}}^{{i}_{N}} \mathrm{s}\mathrm{i}{\mathrm{n}}^{2} \left({\theta }_{jik}\right)\cdot \mathrm{s}\mathrm{i}{\mathrm{n}}^{2} \left({\theta }_{ji{k}^{\prime }}\right)\\ & \quad \cdot {\hat{\beta }}_{ik}^{2}\left({r}_{ik}\right)\cdot {\hat{\beta }}_{i{k}^{\prime }}^{2}\left({r}_{i{k}^{\prime }}\right)\cdot \mathrm{cos}\left({\Delta}{{\Psi}}_{k{k}^{\prime }}\right)\\ & \quad +\frac{1}{2}\sum\limits _{\begin{matrix}\hfill k = {j}_{1}\hfill \\ \hfill k\ne i\hfill \end{matrix}}^{{j}_{N}}\sum\limits _{\begin{matrix}\hfill {k}^{\prime } = k + 1\hfill \\ \hfill {k}^{\prime }\ne i\hfill \end{matrix}}^{{j}_{N}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\left({\theta }_{ijk}\right)\cdot \mathrm{s}\mathrm{i}{\mathrm{n}}^{2} \left({\theta }_{ij{k}^{\prime }}\right)\\ & \quad \cdot {\hat{\beta }}_{jk}^{2}\left({r}_{jk}\right)\cdot {\hat{\beta }}_{j{k}^{\prime }}^{2}\left({r}_{j{k}^{\prime }}\right)\cdot \mathrm{cos}\left({\Delta}{{\Psi}}_{k{k}^{\prime }}\right)\\ & \quad +\frac{1}{2}\sum\limits _{\begin{matrix}\hfill {k}^{\prime }={i}_{1}\hfill \\ \hfill {k}^{\prime }\ne j\hfill \end{matrix}}^{{i}_{N}}\sum\limits _{\begin{matrix}\hfill k={j}_{1}\hfill \\ \hfill k\ne i\hfill \end{matrix}}^{{j}_{N}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\left({\theta }_{ji{k}^{\prime }}\right)\cdot \mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\left({\theta }_{ijk}\right)\\ & \quad \cdot {\hat{\beta }}_{i{k}^{\prime }}^{2}\left({r}_{i{k}^{\prime }}\right)\cdot {\hat{\beta }}_{jk}^{2}\left({r}_{jk}\right)\cdot \mathrm{cos}\left({\Delta}{{\Psi}}_{k{k}^{\prime }}\right),\end{aligned}\end{equation} \tag{ 17 }$

with ${\hat{\beta }}_{ik}^{2}\left({r}_{ik}\right)={p}_{\pi ,i}\cdot {\beta }_{\sigma ,ik}^{2}\left({r}_{ik}\right)-{\beta }_{\pi ,ik}^{2}\left({r}_{ik}\right)$ and $\mathrm{cos}\left({\Delta}{{\Psi}}_{k{k}^{\prime }}\right)=\frac{2{\left(\mathrm{cos}\left({\theta }_{ki{k}^{\prime }}\right)-\mathrm{cos}\left({\theta }_{ji{k}^{\prime }}\right)\cdot \mathrm{cos}\left({\theta }_{jik}\right)\right)}^{2}}{\mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\left({\theta }_{jik}\right)\cdot \mathrm{s}\mathrm{i}{\mathrm{n}}^{2}\left({\theta }_{ji{k}^{\prime }}\right)}-1$ .

In the above analytical functions of BOP, there are altogether four global parameters, seventeen parameters for single element Fe or Bi, and nineteen parameters for the Fe–Bi binary system needed to be determined. The following parameters can be determined according to the principles adopted by Ward et al [20] before the fitting process: the global parameters (ζ₁, ζ₂, ζ₃, ζ₄) and the three-body-dependent parameters g_0,jik for Fe, Bi and Fe–Bi are set as the same values by Ward et al [20] and are shown in tables 1 and 3; the pair-dependent parameters (r₀, r_c, r₁, r_cut) for Fe, Bi and Fe–Bi are displayed in table 2; other selected parameters in table 2 are a_π = c_π = k_σ = p_π = β_π = 0 and f_σ = 0.5 for Fe, and a_π = c_π = 1 for Bi and Fe–Bi.

Table 1. Global and element-dependent BOP parameters for Bi and Fe.

ζ₁	ζ₂	ζ₃	ζ₄	p_π,Bi	p_π,Fe
0.000 01	0.000 01	0.001	0.000 01	0.477 904	0

Table 2. Pair-dependent BOP parameters for Bi, Fe, and Fe–Bi.

	Bi–Bi	Fe–Fe	Fe–Bi
r₀ (Å)	3.1800	2.2452	2.7700
r_c (Å)	3.1800	2.2452	2.7700
r₁ (Å)	3.8200	2.7232	3.3800
r_cut (Å)	4.3000	3.6492	3.9500
n_c	3.497 238	3.289 429	4.441 948
m	1.994 330	1.046 363	1.211 758
n	1.107 474	0.461 114	0.622 392
φ₀ (eV)	1.137 997	1.289 055	0.952 986
β_σ,0 (eV)	0.468 553	1.306 371	1.272 759
β_π,0 (eV)	0.733 087	0	0.620 172
c_σ	0.631 549	0.189 093	1.269 255
f_σ	0.470 068	0.5	0.225 232
k_σ	0.633 972	0	3.904 155
c_π	1	0	1
a_π	1	0	1

The remaining parameters should be obtained from a two-step fitting method as follows [20]. Based on equations (1)–(4), the bond energy (E_b,ij) between atoms i and j can be expressed as:

$\begin{equation}{E}_{\mathrm{b},ij} = {\varnothing }_{ij}\left({r}_{ij}\right) - {\beta }_{\sigma ,ij}\left({r}_{ij}\right)\cdot \frac{2\cdot {\beta }_{\sigma ,0,ij}\cdot {{\Theta}}_{\sigma ,ij}+2\cdot {\beta }_{\pi ,0,ij}\cdot {{\Theta}}_{\pi ,ij}}{{\beta }_{\sigma ,0,ij}}.\end{equation} \tag{ 18 }$

At the equilibrium position, the first derivative ( ${E}_{\mathrm{b},ij}^{\prime }$ ) of E_b,ij equals zero and

$\begin{equation}2\cdot {\beta }_{\sigma ,0,ij}\cdot {{\Theta}}_{\sigma ,ij}+2\cdot {\beta }_{\pi ,0,ij}\cdot {{\Theta}}_{\pi ,ij}={\beta }_{\sigma ,0,ij}\cdot \frac{{\varnothing }_{ij}^{\prime }\left({r}_{ij}\right)}{{\beta }_{\sigma ,ij}^{\prime }\left({r}_{ij}\right)}.\end{equation} \tag{ 19 }$

Substituting equation (19) in equation (18), the bond energy at the equilibrium position can be written as a function of equilibrium bond length:

$\begin{equation}{E}_{\mathrm{b},ij}\left({r}_{ij}\right)={\varnothing }_{ij}\left({r}_{ij}\right)-{\beta }_{\sigma ,ij}\left({r}_{ij}\right)\cdot \frac{{\varnothing }_{ij}^{\prime }\left({r}_{ij}\right)}{{\beta }_{\sigma ,ij}^{\prime }\left({r}_{ij}\right)}.\end{equation} \tag{ 20 }$

Similarly, the second derivative ( ${E}_{\mathrm{b}\cdot ij}^{{\prime\prime}}$ ) of E_b,ij could be derived as follows:

$\begin{equation}{E}_{\mathrm{b}\cdot ij}^{{\prime\prime}}={\varnothing }_{ij}^{{\prime\prime}}\left({r}_{ij}\right)-{\beta }_{\sigma ,ij}^{{\prime\prime}}\left({r}_{ij}\right)\cdot \frac{{\varnothing }_{ij}^{\prime }\left({r}_{ij}\right)}{{\beta }_{\sigma ,ij}^{\prime }\left({r}_{ij}\right)}.\end{equation} \tag{ 21 }$

The atomic volume (V_ij), cohesive energy (E_c,ij), and bulk modulus (B_ij) of a given structure can be then calculated through the following forms: [20]

$\begin{equation}{V}_{ij}=F\cdot {r}_{ij}^{3},\qquad {E}_{\mathrm{c},ij}=\frac{Z\cdot {E}_{\mathrm{b},ij}}{2},\qquad {B}_{ij}\left({r}_{ij}\right)=\frac{Z\cdot {E}_{\mathrm{b}\cdot ij}^{{\prime\prime}}}{18F\cdot {r}_{ij}},\end{equation} \tag{ 22 }$

where F is a structural volume factor, and Z is the atomic coordination.

For the first-step fitting, equations (20) and (22) are used to determine the pairwise GSP parameters (∅₀, m, n, and n_c) through the fitting of E_c and B of each selected structure at the equilibrium volume from experiments or first principles calculations. Specifically, the structures of A7, BCC, FCC, HCP, and SC (simple cubic) are selected for Bi, the BCC, FCC, HCP, and SC structures are chosen for Fe, and the B1, B2, B3, B4, and C3 structures for Fe–Bi.

The second step of fitting is to obtain the remaining parameters (β_σ,0, β_π,0, c_σ, f_σ, k_σ, b_σ, u_σ), i.e., the derived parameters (∅₀, m, n, and n_c) in the first step are substituted into the right-hand side of equation (19) and treated as new reference values for the left-hand side of equation (19) during the fitting process. In addition, equations (18) and (22) are also used in the second step fitting, and the input data for equations (18) and (22) are the E_c values of the following structures from experiments or first principles calculations: A7 and BCC for Bi, BCC for Fe, and B1 and B2 for Fe–Bi.

It should be noted that the present first-principles calculations are performed through the Vienna ab initio simulation package (VASP) software [32]. The interaction between ion and valence electrons is described by the projector-augmented wave method [33, 34], and the exchange and correlation effect is expressed by the generalized gradient approximation [35–37]. The plane-wave cutoff energy is 500 eV for all the calculations. In addition, the temperature-smearing method of Methfessel–Paxton [38] and the tetrahedron method with Blöchl corrections [39] are performed for relaxation and static calculations during k space integration, respectively. The unit cells of A7, BCC, FCC, HCP, SC, B1, B2, B3, B4, and C3 structures contain 6, 2, 8, 2, 1, 8, 2, 8, 4, and 6 atoms, respectively, and their k-points are 15 × 15 × 7, 11 × 11 × 11, 11 × 11 × 5, 11 × 11 × 7, 11 × 11 × 11, 11 × 11 × 11, 11 × 11 × 11, 11 × 11 × 11, 11 × 11 × 7, and 11 × 11 × 11, respectively.

During the process of fitting, the properties of Bi with the A7 structure as well Fe with BCC and FCC structures have the higher weights than others. After a series of fitting, tables 1–3 summarize the optimized parameters of Bi, Fe, and Fe–Bi, and table 4 lists the obtained a, c, E_c, and B values of Bi, Fe, and FeBi with the selected structures from the present BOP, first principles calculations, and experiments [40–45]. It should be mentioned that the values of cohesive energy E_c from first principles calculations were shifted to experimental values according to the calculated energies differences for the sake of comparation. One can observe from table 4 that the constructed BOP has a generally nice performance in reproducing the E_c, and B values of Bi, Fe, and FeBi from experiments [40–45] and first principles calculation. The validation of the constructed Fe–Bi potential will be shown later in section 4.

Table 3. Three-body-dependent BOP parameters for Bi, Fe, and Fe–Bi.

	BiBiBi	FeFeFe	BiBiFe	BiFeBi	BiFeFe	FeBiFe
g_0,jik	1	1	1	1	1	1
b_σ,jik	0.488 041	0.920 095	0.862 671	0.557 518	0.698 533	0.830 309
u_σ,jik	−0.114 961	−0.074 563	−0.121 116	−0.069 474	−0.189 584	−0.179 300

Table 4. Lattice parameters (a, c), cohesive energies (E_c), and bulk moduli (B) of Bi, Fe and FeBi with several selected structures from the present potential (BOP), first principles calculations (VASP), and experiments [40–45].

Phase	Structure	a, c (Å)			E_c (eV/atom)			B₀ (GPa)
Phase	Structure	BOP	VASP	Exp.	BOP	VASP	Exp.	BOP	VASP	Exp.
Bi	A7	4.563	4.540	4.546 [41]	−2.1803	−2.180	−2.18 [42]	55	35	31 [40]
	A7	11.827	11.577	11.865 [41]	−2.1803	−2.180	−2.18 [42]	55	35	31 [40]
	SC	3.260	3.219		−2.1798	−2.144		58	50
	FCC	4.988	4.949		−1.8631	−2.059		87	54
	BCC	3.960	3.912		−1.9328	−2.085		88	54
	HCP	3.520	3.547		−1.8613	−2.083		87	52
	HCP	5.780	5.792		−1.8613	−2.083		87	52
Fe	BCC	2.866	2.832	2.866 [45]	−4.2832	−4.280	−4.28 [42]	163	185	169 [43]
	FCC	3.680	3.539	3.640 [45]	−4.2750	−4.179		133	172	133 [44]
	SC	2.160	2.373		−3.4324	−3.515		120	103
	HCP	2.602	2.430		−4.2762	−4.189		135	290
	HCP	4.249	3.968		−4.2762	−4.189		135	290
FeBi	B1	5.737	5.594		−2.4868	−2.562		65	62
	B2	3.431	3.415		−2.4813	−2.640		76	75
	B3	6.151	6.161		−2.3085	−2.299		38	36
	B4	4.351	4.353		−2.3085	−2.297		34	38
	B4	7.106	7.108		−2.3085	−2.297		34	38
Fe₂Bi	C3	6.100	5.982		−1.6124	−1.757		12	32
FeBi₂	C3	6.019	5.636		−2.0451	−2.156		32	48

3. Simulation methods

Based on the newly constructed BOP of Fe–Bi, the present MD simulations are performed by means of the well-known software of large-scale atomic/molecular massively parallel simulator package [46]. Accordingly, the Fe_100−xBi_x solid solutions are constructed within the entire composition range of 0 ⩽ x ⩽ 100 with an interval of x = 2. Two structures of BCC and A7 are selected at each composition of Fe_100−xBi_x, and the unit cells of 68 × 68 × 68 Å³ (16 000 atoms) and 81.6 × 86.4 × 82.7 Å³ (16 632 atoms) are chosen for BCC and A7 structures, respectively. The x, y, and z directions of the unit cell are set along [100], [010], and [001] for the BCC structure, respectively, while [ $2\bar{1}\bar{1}0$ ], [ $01\bar{1}0$ ], and [0001] for the A7 structure, respectively. For each composition of Fe_100−xBi_x, the unit cell is formed by randomly substituting Fe (or Bi) atoms by Bi (or Fe) atoms. Periodic boundary conditions are added at three directions of the unit cell, and each solid solution is relaxed at 0.1 K under the isothermal-isobaric (NPT) ensemble for 50 ps and optimized through the conjugate gradient method. It should be noted that for BCC or A7 structure at each composition of Fe_100−xBi_x, five samples are considered for random substitutions, and the data shown in the following tables and figures are average values. The relative errors of the calculated lattice parameters, cohesive energies, elastic constants and moduli, surface energies, heat capacities, and heat of formation are within 0.001%, 10⁻⁶ %, 2%, 1%, 1%, and 5%, respectively.

4. Result and disscusion

It is of importance to validate the newly constructed BOP of Bi–Bi, Fe–Fe, and Fe–Bi. After a series of MD simulations, it will show that the present BOP of Fe–Bi could correctly predict the ground state A7 structure of Bi, reveal the structural stability of Fe–Bi phases, and find out the effects of Bi on mechanical properties of BCC Fe. The simulated results from the present potential are in good agreement with experimental data in the literature, and the present BOP of Fe–Bi has a better performance to reflect the intrinsic properties of the binary Fe–Bi system than the two LJ potentials and an EAM potential have been constructed for Bi and Fe–Bi [15–18].

4.1. Ground-state structure and properties of Bi

It is well known that the ground-state structure of Bi is a rhombohedral structure of A7 (space group No. 166, R $\bar{3}$ m, Pearson symbol hR2) with a special rhombohedral angle (θ_rh) of 57.23° [41]. This structure can be also recognized as a hexagonal structure with a unique internal parameter of u = 0.2340 as a result of its layered atomic arrangement [47]. The A7 structure of Bi and its equivalent hexagonal structure are shown vividly in figure 1. Such a special angle of θ_rh and an internal parameter of u make it difficult to accurately reproduce the ground-state structure of Bi by an n-body potential. Probably due to this difficulty, the A7 structure has been simply disregarded by the EAM and two LJ potentials of Bi in the literature and the BCC structure of Bi is used for the fitting instead [15–18]. In the present study, BOP has been intentionally selected to construct the potential of Bi due to its suitability for directional bonding, and it is therefore of importance to see whether or not the constructed BOP of Bi could reflect the ground-state structure and properties of Bi. The properties from the EAM and two LJ potentials of Bi already in the literature are also included for the sake of comparison.

**Figure 1.** Rhombohedral (A7) and hexagonal unit cell of Bi bulk. The brown and green atoms distinguish the atoms in the two layers of each bilayer.
Download figure:
Standard image High-resolution image

Consequently, figure 2 shows the comparison of energy curves of Bi with the structures of A7, BCC, FCC, HCP, and SC obtained from the potentials of the present BOP, EAM [17, 18], as well as LJ of Maulana [16] and Liu [15]. It can be seen clearly from figure 2(a) that the A7 structure is predicted by the present BOP to possess the lowest energy among the selected structures, and is therefore regarded as the ground-state structure of Bi, which agrees well with the ground-state A7 structure of Bi from experimental observations in the literature [41, 42]. On the contrary, the energy curves in figure 2 from the EAM and LJ potentials of Bi reveal that the HCP or FCC structure has the lowest energy, which seems contradictory to the ground-state A7 structure by experiments [41, 42].

Table 5 summarizes the derived cohesive energy (E_c), lattice parameters (a, a_rh and c/a), internal parameter (u), and rhombohedral angle (θ_rh) of the A7 structure of Bi obtained from the present BOP, experiments [41, 42], EAM potentials [17, 18], and LJ potentials of Maulana [16] and Liu [15]. One could discern from this table that the present BOP could accurately predict the structural parameters of the ground-state A7 structure of Bi. For instance, the a_rh, u, and θ_rh of A7 from the present study are 4.741 Å, 0.2369, and 57.53°, respectively, which match well with the corresponding experimental values of 4.746 Å, 0.2340, and 57.23°. Unfortunately, the EAM and LJ potentials of Bi [15–18] could not correctly reproduce the cohesive energy and structural parameters of the A7 structure due to the big differences of the calculated and experimental values shown in table 5.

Table 5. Cohesive energy (E_c), lattice parameters (a, a_rh and c/a), internal parameter (u), and rhombohedral angle (θ_rh) of the A7 structure of Bi obtained from the present potential (BOP), experiments [41, 42], EAM potentials [17, 18], and the LJ potentials of Maulana [16] and Liu [15].

	E_c (eV/atom)	a (Å)	a_rh (Å)	u	c/a	θ_rh (°)
BOP	−2.1803	4.563	4.741	0.2369	2.592	57.53
Exp. [41, 42]	−2.18	4.546	4.746	0.2340	2.610	57.23
EAM [17, 18]	−2.0693	5.618	4.138	0.2491	1.229	86.79
LJ [16]	−0.6123	4.889	3.457	0.2500	1.225	90
LJ [15]	−0.4803	5.326	3.767	0.2500	1.225	90

In addition, the elastic constants of the A7 structure of Bi are calculated through the relationship of energy and strain [48], and the obtained results are displayed in table 6. The C₁₁, C₁₄, and C₄₄ values of the A7 structure from the present BOP are in accordance with experimental measurements [40], while the present potential overestimates the C₁₂, C₁₃, and C₃₃ values. One could also observe that the elastic constants from the EAM potential [17, 18] are much bigger than the experimental values [40], and that the present BOP has a much nicer agreement with experiments than the EAM potential. It should be pointed out that the A7 structure of Bi from the two LJ potentials [15, 16] becomes mechanically unstable and its elastic constants are thus unavailable.

Table 6. Elastic constants C_ij (GPa) and bulk moduli B (GPa) of A7 structure of Bi as well as BCC and FCC structures of Fe from the present potential (BOP), the present first principles calculation (VASP), experiments [40, 43, 44], EAM potentials of Bi [17, 18] and Fe [51], and the LJ potentials of Maulana [16] and Liu [15].

Phase	Structure	C_ij	BOP	VASP	Exp. [40, 43, 44]	EAM [17, 18, 51]	LJ [16]	LJ [15]
Bi	A7	C₁₁	66	84	69.3	106	Unstable	Unstable
		C₁₂	47	24	24.5	53	Unstable	Unstable
		C₁₃	52	35	25.4	58	Unstable	Unstable
		C₁₄	4	21	8.4	−28	Unstable	Unstable
		C₃₃	73	58	40.4	168	Unstable	Unstable
		C₄₄	16	41	13.5	195	Unstable	Unstable
Fe	BCC	C₁₁	223	258	226	229	1187	679
		C₁₂	133	149	140	135	499	387
		C₄₄	100	96	116	116	499	146
Fe	FCC	C₁₁	196	253	154	107	876	679
		C₁₂	101	131	122	98	499	387
		C₄₄	98	106	77	80	499	387

The surface energy and heat capacity of the A7 structure of Bi are also calculated by MD simulation with the present BOP, and the results are shown in table 7. The corresponding experimental value [49] as well as available calculated data from first principles calculation [50] and EAM potential [17, 18] are also listed in table 7 for comparison. The present surface energy (0.17 J m⁻²) of the (111) surface of Bi matches well with the corresponding value of 0.16 J m⁻² from first principles calculations [50], whereas such a value of 0.44 J m⁻² from the EAM potential [17, 18] seems much bigger. Moreover, the volumetric heat capacity of the A7 structure of Bi at 283 K from the present BOP is compatible with the corresponding values from experiments [49] and the EAM potential [17, 18].

Table 7. Surface energy E_s (J m⁻²) and heat capacity C_v (J·mol⁻¹·K⁻¹) at 283 K of A7 Bi and BCC Fe from the present potential (BOP), experiments [49, 52, 53], first principles calculation [50, 54], EAM potentials of Bi and Fe [17, 18, 51], and the LJ potentials of Maulana [16] and Liu [15]. The experimental E_s of Fe is surface energy of polycrystalline sample.

	Bi		Fe
	E_s (111)	C_v	E_s (100)	E_s (110)	E_s (111)	C_v
BOP	0.17	28.5	2.06	1.89	2.34	25.3
Exp. [49, 52, 53]	—	25.9	2.36			25.6
Cal. [50, 54]	0.16	—	2.29	2.27	2.52	—
EAM [17, 18, 51]	0.44	28.2	1.70	1.43	1.94	25.1
LJ [16]	—	—	5.89	5.65	6.28	24.2
LJ [15]	—	—	4.68	4.48	4.98	24.3

From the above simulated results, one can see clearly that the constructed BOP of Bi is not only able to successfully predict the ground-state A7 structure of Bi and its structural parameters, but also capable of reproducing the elastic constants, surface energy, and heat capacity of the A7 structure of Bi. These nice agreements with experiments validate the relevance of the present BOP of Bi. The EAM and LJ potentials of Bi already in the literature [15–18], however, are unable to reflect the ground-state A7 structure of Bi and its structural parameters. It should be noted that the present BOP potential of Bi has been constructed, and this BOP potential of Bi has been therefore used to derive the Fe–Bi cross potential, the validation of which will be presented in the following sections.

Additionally, it is of interest to see the performance of the BOP of Fe. After a series of calculations, table 6 lists the obtained elastic constants (C₁₁, C₁₂, C₄₄) of Fe with the BCC and FCC structures, and table 7 shows the calculated surface energies (E_s), and heat capacity (C_v) of BCC Fe. One can discern clearly from these two tables that the C₁₁, C₁₂, C₄₄, E_s, and C_v values of Fe from the present BOP are compatible with corresponding experimental and calculated data in the literature [43, 44], suggesting that the constructed BOP of Fe is realistic to reflect the intrinsic feature of Fe. As shown also in table 7, the EAM potential [51] has a similar performance to reproduce the elastic constants and heat capacity of Fe as the present BOP, while underestimates the surface energy of Fe. On the other hand, the two LJ potentials [15, 16] considerably overestimate the elastic constants and surface energy of Fe.

4.2. Structural stability of Fe–Bi phases

As related before, Fe–Bi is an immiscible binary system with a positive heat of formation (+18 kJ mol⁻¹) [1], while metastable Fe–Bi phases with the BCC and A7 structures have been fabricated by ion beam mixing [8], magnetron sputtering [3, 4, 8], electron-beam evaporation [2, 9, 10], and mechanical alloying [11], etc. Interestingly, amorphous Fe–Bi phases have been also discovered experimentally when the Bi composition is in the range of 40–60 at.% and 10–30 at.% [8–10]. As to theoretical studies about the metastable Fe–Bi phases, however, there is not any investigation so far in the literature. Based on the newly constructed BOP of Fe–Bi, the present study is therefore dedicated to find out the structural stability of the binary Fe–Bi system.

Accordingly, the unit cells of Fe_100−xBi_x solid solutions with the BCC and A7 structures are constructed, respectively, within the entire composition range at an x interval of 2. Each Fe_100−xBi_x solid solution is obtained through a random substitution of Fe (or Bi) atoms with Bi (or Fe). After relaxation, it is found that the Fe_100−xBi_x solid solutions could keep the BCC or A7 structure when x is within a certain range. On the other hand, the BCC or A7 structure of Fe_100−xBi_x solid solutions would become unstable, and will completely or partially transform to amorphous structure when x exceeds a specific composition. The revealed composition ranges of the BCC, A7, and amorphous structures of Fe_100−xBi_x will be described in the following paragraphs.

The heat of formation (ΔH_f) of each Fe_100−xBi_x phase with the BCC, A7, or amorphous structure is calculated by the following form:

$\begin{equation}{\Delta}{H}_{\mathrm{f}}={E}_{\mathrm{F}{\mathrm{e}}_{100-x}\mathrm{B}{\mathrm{i}}_{x}}-x{E}_{\mathrm{B}\mathrm{i}}-(100-x){E}_{\mathrm{F}\mathrm{e}},\end{equation} \tag{ 23 }$

where ${E}_{\mathrm{F}{\mathrm{e}}_{100-x}\mathrm{B}{\mathrm{i}}_{x}}$ represents the total energy of each Fe_100−xBi_x phase, E_Bi and E_Fe are totals energies of the ground-state A7 and BCC structures for Bi and Fe, respectively. It should be mentioned that the term of the product of pressure and volume is not included in equation (23) due to the zero pressure maintained during the present MD simulations. After the calculations, figure 3 shows the heats of formation of the Fe_100−xBi_x phases obtained from the present study, the LJ potentials of Maulana [16] and Liu [15], as well as the thermodynamic model of Miedema et al [1]. In addition, figures 4 and 5 display, as typical examples, the atomic configurations and radial distribution functions of several Fe_100−xBi_x phases, respectively, and figure 6 shows the partial radial distribution functions of Fe–Fe, Fe–Bi, and Bi–Bi in the amorphous structure of Fe₅₀Bi₅₀.

**Figure 3.** Heats of formation of Fe_100−xBi_x phases obtained from the present study, the LJ potentials of Maulana [16] and Liu [15], as well as the thermodynamic model of Miedema *et al* [1].
Download figure:
Standard image High-resolution image

**Figure 4.** Atomic configurations of (a) BCC + amorphous of Fe₉₀Bi₁₀, (b) amorphous of Fe₅₀Bi₅₀, (c) A7 + amorphous of Fe₃₀Bi₇₀, and (d) A7 of Fe₁₄Bi₈₆. The blue and brown balls represent Fe and Bi atoms, respectively. The x, y, and z directions are along [100], [010], [001] of the BCC structure, and along [ $2\bar{1}\bar{1}0$ ], [ $01\bar{1}0$ ], [0001] of the A7 structure, respectively.
Download figure:
Standard image High-resolution image

**Figure 4.** Atomic configurations of (a) BCC + amorphous of Fe₉₀Bi₁₀, (b) amorphous of Fe₅₀Bi₅₀, (c) A7 + amorphous of Fe₃₀Bi₇₀, and (d) A7 of Fe₁₄Bi₈₆. The blue and brown balls represent Fe and Bi atoms, respectively. The x, y, and z directions are along [100], [010], [001] of the BCC structure, and along [ $2\bar{1}\bar{1}0$ ], [ $01\bar{1}0$ ], [0001] of the A7 structure, respectively.
Download figure:
Standard image High-resolution image

**Figure 5.** Radial distribution functions of (a) BCC Fe, (b) Fe₉₀Bi₁₀, (c) Fe₅₀Bi₅₀, (d) Fe₃₀Bi₇₀, and (e) A7 Bi.
Download figure:
Standard image High-resolution image

**Figure 6.** Partial radial distribution functions of Fe–Fe, Fe–Bi, and Bi–Bi in the amorphous structure of Fe₅₀Bi₅₀.
Download figure:
Standard image High-resolution image

Some features of structural stability of Fe_100−xBi_x phases could be discerned from figures 3–6. First of all, the heats of formation of the energetically favorable Fe_100−xBi_x phases obtained from the present study are compatible with those from the thermodynamic model of Miedema et al [1], although some apparent discrepancy exists between the present and thermodynamic values within some composition ranges shown in figure 3. This discrepancy would be probably due to the fact that crystal structure was not considered during thermodynamic calculations [1]. On the other hand, the two LJ potentials of Fe–Bi [15, 16] significantly overestimate the heats of formation of BCC Fe_100−xBi_x phases, and could not simulate the A7 structure of Fe_100−xBi_x due to its instability.

Secondly, the Fe_100−xBi_x phases could keep the BCC and A7 structures when 0 ⩽ x < 10 and 86 ⩽ x ⩽ 100, respectively. By means of mechanical alloying, Grigorieva et al discovered that Fe₉₇Bi₃ has a single BCC structure [11], and the composition of Fe₉₇Bi₃ is just located within the predicted composition range (0 ⩽ x < 10) of BCC Fe_100−xBi_x phases. Moreover, Bi atoms are experimentally found to dissolve in the BCC lattice of Fe, and Fe atoms are discovered to substitute the Bi atoms in the A7 structure [3, 9], which are also in accordance with the above simulated results.

Thirdly, the BCC Fe_100−xBi_x phases become partially amorphous when 10 ⩽ x < 26, i.e., BCC + amorphous. As a typical example to vividly show the mixture of BCC and amorphous, one can see from figure 4(a) that some disordered regions called amorphous have been formed in the atomic configurations of Fe₉₀Bi₁₀, while its BCC structure can be still apparently discerned. In addition, the radial distribution functions of Fe₉₀Bi₁₀ displayed in figure 5(b) also have typical features of both BCC and amorphous structures, i.e., the original sharp peaks of BCC shown in figure 5(a) have been broadened or disappeared when the Bi composition reaches the critical point of 10 at.%.

Fourthly, on the Bi-rich side, partial amorphization (A7 + amorphous) of the A7 Fe_100−xBi_x phases could be also realized when the Bi composition is within the range of 70 ⩽ x < 86. Compared with the A7 structure of Fe₁₄Bi₈₆ in figure 4(d), the atomic configuration of Fe₃₀Bi₇₀ in figure 4(c) exhibits a largely disordered pattern, whereas the A7 structure could be observed to some extent. Such a mixture of A7 + amorphous of Fe₃₀Bi₇₀ could also be seen in terms of pair correlation function in figure 5(d). Interestingly, the partial amorphization (A7 + amorphous) of FeBi phases was indeed observed experimentally [2, 9], which is in nice agreement with the above simulated results from the present BOP of Fe–Bi.

Fifthly, a uniform amorphous structure of Fe_100−xBi_x could be formed when the composition of Bi is located in the range of 26 ⩽ x < 70. Figure 4(b), as a typical example, displays the atomic configurations of the obtained amorphous structure of Fe₅₀Bi₅₀. A completely disordered pattern has been formed in Fe₅₀Bi₅₀, and the BCC or A7 structure could no longer be discerned in figure 4(b). Moreover, figure 5(c) also shows the pair correlation function of Fe₅₀Bi₅₀, which has a typical feature of amorphous with only a short-range order.

It is of importance to compare the above simulation results regarding amorphous Fe_100−xBi_x with similar experimental observations in the literature. Chen et al discovered that the Fe_100−xBi_x phase becomes a uniform amorphous structure when Bi is in the range of 40–60 at.% and 45–60 at.% after ion-beam mixing [8] and electron beam codeposition [9], respectively. Such experimentally-observed composition ranges are just located within the above predicted range of 26 ⩽ x < 70 for amorphous structure of Fe_100−xBi_x, and it seems understandable that the present theoretical composition range of amorphous formation seems a little bigger. By means of experiments, Forester et al also identified the partial amorphization of Fe_100−xBi_x when 10 ⩽ x ⩽ 30 [10], which is in good agreement with the present predicted composition of 10 ⩽ x < 26 for the BCC + amorphous structure of Fe_100−xBi_x. Furthermore, the atomic distance of Fe–Bi in the amorphous structure of Fe₅₀Bi₅₀ is experimentally measured to be 2.8 Å [8], which matches well with the present simulated value of 2.8 Å shown in figure 6.

All the above points indicate that the present MD simulation based on the newly constructed BOP of Fe–Bi is able to reveal the structural stability of the binary Fe–Bi system in terms of metastable phase formation. The predicted BCC, BCC + amorphous, amorphous, A7 + amorphous, and A7 structures of Fe_100−xBi_x phases within the entire composition range agree well with similar experimental observations in the literature [2, 3, 8–11]. These nice agreements therefore provide a strong support for the validation of the newly constructed BOP of Fe–Bi. Nevertheless, it should be pointed out that the two LJ potentials [15, 16] are unable to reflect the structural stability of various Fe–Bi phases.

4.3. Effects of Bi on mechanical properties of BCC Fe

It is commonly believed that Bi and its alloys are promising candidates of spallation target and coolant with superior properties in nuclear power systems. The interaction between Fe and Bi, however, has been a challenging issue among scientists and engineers, as Bi and its alloys are found to corrode the structural materials, and bring about the so-called LME of ferrite/martensite stainless steels, the mechanism of which has been still controversial during the past decades [12–14]. In the present MD simulations, it is therefore of importance to find out the effect of Bi on elastic properties and stress–strain relationships of BCC Fe based on the newly constructed Fe–Bi potential.

After a series of calculations, figure 7 shows the obtained elastic constants (C₁₁, C₁₂, C₄₄), elastic moduli (B, G, E), and G/B ratios of BCC Fe_100−xBi_x phases as a function of Bi composition. It should be noted that five compositions of Fe_100−xBi_x with an interval of x = 2 have been selected in figure 7 and the maximum composition of Bi is 8 at.%, as the stable range of the BCC structure of Fe_100−xBi_x is 0 ⩽ x < 10 as derived before. One can observe from figure 7 that the addition of Bi has an important effect to reduce the elastic constants and elastic moduli of BCC Fe, and such a decrease become more apparent with the increase of Bi composition. This effect of Bi would be probably ascribed to the much lower elastic constants of Bi than those of Fe as shown in table 6.

The G/B ratios of these five BCC Fe_100−xBi_x phases are also calculated and are displayed in figure 7. It should be pointed out that the G/B ratio has been universally regarded as a criterion to express the brittle/ductile behavior of a materials, i.e., a smaller G/B less than the critical value of 0.57 signifies a more ductile feature, and vice versa. Interestingly, the G/B ratios shown in figure 7 decrease with the increase of the Bi composition, suggesting that the addition of Bi could make it more ductile for BCC Fe, which is apparently contradictory to the experimental observation of the brittle behavior of BCC steels due to the effect of Bi [12–14]. Such a disagreement indicates that the G/B ratio in the elastic range should be unsuitable to describe the effect of Bi on the brittle/ductile behavior of BCC Fe, and that plastic deformation would probably occur during the experimental brittle process of the Fe–Bi system [12–14].

Furthermore, figure 8 displays, as typical examples, the stress–strain curves of BCC Fe_100−xBi_x phases with four Bi compositions (x = 0.0065, 0.1, 1, and 2) under a uniaxial tensile loading along the [111] direction. Each tensile loading is performed at 0.1 K with a constant engineering strain rate of 1 × 10⁹ s⁻¹ and a timestep of 0.001 ps under the NPT ensemble. The corresponding stress–strain curve of BCC Fe is also included in figure 8 for the sake of comparison. It should be pointed out that the maximum Bi composition shown in figure 8 is 2 at.%, since the BCC structure of Fe_100−xBi_x is found to become unstable with an incorrect stress–strain curve under a uniaxial loading when x is bigger than 2.

One can discern clearly from figure 8 that the stress–strain curve of BCC Fe_100−xBi_x phase when x = 0.0065 is apparently lower with a much smaller critical strain than that of pure BCC Fe, implying that a very small percentage (0.0065 at.%) of Bi could considerably reduce the tensile strength and ductility of BCC Fe. Moreover, both tensile strength and critical strain of BCC Fe_100−xBi_x phases decrease with the increase of the Bi composition. In other words, the addition of Bi would significantly increase the brittleness of BCC Fe, which is in good agreement with similar experimental observations in the literature [12–14]. Such an agreement would provide further evidence for the relevance of the newly constructed BOP of Fe–Bi.

It is of interest to find out the fundamental reason why a very small percentage (0.0065 at.%) of Bi could considerably reduce the ductility of BCC Fe. Accordingly, figure 9 displays the atomic configuration and dislocation lines of the Fe_1−xBi_x phase (x = 0.0065 at.%) under uniaxial tensile loading along the [111] direction at the strain of 16.6%. As shown in figure 9, the addition of Bi has caused a considerable local distortion of the nearby Fe atoms and some dislocations lines of 1/2 ⟨111⟩ around the Bi atoms at the strain of 16.6%. Pure BCC Fe, however, only possesses a uniform deformation without any local distortion and dislocation lines during the uniaxial tensile loading (figures not shown). Such a local distortion and dislocations lines around the Bi atom would be probably attributed to the bigger atomic radius of Bi and the repulsive interaction between Bi and Fe, which could intrinsically bring about the much lower tensile strength and ductility of BCC Fe as shown in figure 8.

**Figure 9.** (a) Atomic configuration and (b) (c) dislocation lines of the BCC Fe_100−xBi_x phase (x = 0.0065 at.%) under uniaxial tensile loading along the [111] direction at the strain of 16.6%. The brown atom refers to the Bi atom and the figures are visualized with the software of OVITO. The x, y, and z directions are along [ $0\bar{1}1$ ], [ $2\bar{1}\bar{1}$ ], and [111], respectively.
Download figure:
Standard image High-resolution image

**Figure 9.** (a) Atomic configuration and (b) (c) dislocation lines of the BCC Fe_100−xBi_x phase (x = 0.0065 at.%) under uniaxial tensile loading along the [111] direction at the strain of 16.6%. The brown atom refers to the Bi atom and the figures are visualized with the software of OVITO. The x, y, and z directions are along [ $0\bar{1}1$ ], [ $2\bar{1}\bar{1}$ ], and [111], respectively.
Download figure:
Standard image High-resolution image

5. Conclusions

In summary, a new Fe–Bi potential of BOP has been constructed and has been proven to be more realistic than the LJ potentials of Fe–Bi already in the literature. Based on this BOP, the present MD simulation has been used to accurately reproduce the ground-state structure (A7) of Bi as well as its structural parameters and properties. Simulations also predict that the structures of BCC, BCC + amorphous, amorphous, A7 + amorphous, and A7 of Fe_100−xBi_x are energetically more favorable with less heat of formation when 0 ⩽ x < 10, 10 ⩽ x < 26, 26 ⩽ x < 70, 70 ⩽ x < 86, and 86 ⩽ x ⩽ 100, respectively. In addition, it is found that the addition of a very small percentage of Bi should have an important effect to significantly decrease the tensile strength of BCC Fe and considerably increase its brittleness. The obtained results are discussed and compared with experimental or theoretical results in the literature and the agreements between them are fairly nice.

Data availability statement

All data that support the findings of this study are included within the article (and any supplementary files).

Atomistic modelling of the immiscible Fe–Bi system from a constructed bond order potential

Article metrics

Submit

Permissions

Author e-mails

Author affiliations

Author notes

ORCID iDs

Dates

Peer review information

Abstract

1. Introduction

2. Construction of Fe–Bi bond-order potential

3. Simulation methods

4. Result and disscusion

4.1. Ground-state structure and properties of Bi

4.2. Structural stability of Fe–Bi phases

4.3. Effects of Bi on mechanical properties of BCC Fe

5. Conclusions

Data availability statement

Atomistic modelling of the immiscible Fe–Bi system from a constructed bond order potential

Article metrics

Submit

Permissions

Share this article

Author e-mails

Author affiliations

Author notes

ORCID iDs

Dates

Peer review information

Abstract

1. Introduction

2. Construction of Fe–Bi bond-order potential

3. Simulation methods

4. Result and disscusion

4.1. Ground-state structure and properties of Bi

4.2. Structural stability of Fe–Bi phases

4.3. Effects of Bi on mechanical properties of BCC Fe

5. Conclusions

Data availability statement