A Modified Embedded-Atom Method Interatomic Potential for the Mg–Mn Binary System

Abstract

Magnesium, the lightest structural metal, needs to improve its strength and formability at room temperature for commercialization. Manganese is an alloying element that can improve the strength of magnesium through forming intermetallic compounds. However, its effect on the formability of magnesium has not yet been clarified. To identify the effect, an interatomic potential for the Mg–Mn binary system has been developed on the basis of the second nearest-neighbor modified embedded-atom method formalism. The Mg–Mn potential reproduces the structural, elastic, and thermodynamic properties of the compound and solution phases of its associated alloy system, consistent with experimental data and higher-level calculations. The applicability of the developed potential is demonstrated by calculating the generalized stacking fault energy for various slip systems of the Mg–Mn alloy.

Fig. 4 Calculated GSFE (mJ/m²) curves of pure Mg and Mg–Mn, models on the (a) basal plane along the $\langle 10\bar{1}0\rangle $ direction, (b) basal plane along the $\langle 11\bar{2}0\rangle $ direction, (c) prismatic plane along the $\langle 11\bar{2}0\rangle $ direction, (d) pyramidal I plane along the $\langle 11\bar{2}0\rangle $ direction, (e) pyramidal I plane along the $\langle 11\bar{2}3\rangle $ direction, and (f) pyramidal II plane along the $\langle 11\bar{2}3\rangle $ direction, using the present 2NN MEAM potential.

1. Introduction

Light-weight, high-strength, and high-formability at room temperature (RT) are the main keywords in recent material development, especially in the automotive industry.^1,2) Magnesium (Mg), the lightest structural metal, satisfies the first keyword well.³⁾ Despite this, Mg has lower strength and formability at RT than steels or aluminum (Al) alloys.⁴⁾ Hence, both these properties need to be improved for its commercialization.

A general method to improve the strength of Mg is precipitation hardening.⁵⁾ Precipitates can be formed by adding alloying elements that form intermetallic compounds with Mg or other frequently added alloying elements, such as Al and zinc (Zn). Among the alloying elements frequently added to Mg, manganese (Mn) satisfies the above conditions. Mn does not form any intermetallic compounds with Mg,⁶⁾ but forms seven and five intermetallic compounds with Al⁷⁾ and Zn,⁸⁾ respectively. These Mn-containing intermetallic compounds of various compositions and structures could improve the strength of Mg by forming many precipitates.

Despite their strengthening abilities, the effect of Mn on the RT formability of Mg is still unclear. It has been reported that this effect of Mn varies depending on the deformation process; an extruded Mg–1 mass% Mn alloy exhibits higher ductility and formability at RT than pure Mg,⁹⁾ but a hot rolled Mg–1 mass% Mn alloy shows lower ductility at RT than pure Mg.¹⁰⁾ Because of the different experimental results, a more fundamental study is necessary to reveal the effect of Mn on Mg formability.

Previous studies have reported that the Mg formability depends on the slip, twin, and recrystallization behaviors.⁵⁾ These three behaviors are fundamentally related to dislocations that cause a slip, contribute to nucleation and growth of twins,^11,12) and nucleate subgrains during recrystallization.¹³⁾ Hence, we need to investigate the effect of Mn on the dislocation behavior of Mg to further understand its fundamental effect on the Mg formability. The dislocation behavior is a dynamic material phenomenon at an atomic scale that changes with time. Therefore, its observation can be performed using tools that allow a continuous nanoscale analysis of the samples, such as an in situ transmission electron microscope (TEM), a first-principles calculation, and a large-scale atomistic simulation. The in situ TEM can observe the motion of the dislocations in a sample, but it is time- and cost-consuming to observe all kinds of Mg dislocations, such as basal, prismatic, pyramidal I, and pyramidal II dislocations. The first-principles calculation can analyze a basal dislocation with a simple core structure, but analyzing non-basal dislocations with complex core structures is a challenge. The large-scale atomistic simulation is comparatively more appropriate for investigating those various Mg dislocations because it can analyze the movement of millions of atoms at an atomic scale. For this reason, the large-scale atomistic simulation method was chosen to investigate the dislocation behavior of the Mg–Mn alloys.

The large-scale atomistic simulation can be conducted based on (semi-)empirical interatomic potentials describing relevant alloy systems. This means that an interatomic potential for the Mg–Mn binary system is required to perform the dislocation investigation. However, this potential has not yet been developed. Therefore, we should develop the interatomic potential describing the Mg–Mn binary system.

This study aims to develop the interatomic potential for the Mg–Mn binary system, as a prior study to investigate the dislocation behavior of the Mg–Mn alloy using the large-scale atomistic simulation, and finally reveal the effect of Mn on the RT formability of Mg. This paper begins by introducing the procedure to determine the potential parameters. Next, we confirm the reliability of the developed potential by computing fundamental material properties and comparing them with the data available from the literature.

2. Determination of Interatomic Potential Parameters

Developing an interatomic potential means determining the model parameter values that compose the interatomic potential formalism for its associated system. The formalism used in this study is the second nearest-neighbor modified embedded-atom method (2NN MEAM) formalism,¹⁴⁾ which has been widely used to develop interatomic potentials for Mg alloy systems.^15–18) Under the 2NN MEAM formalism, a potential for a binary system comprises 13 binary model parameters, and 28 unary model parameters that compose the potentials for the two constituent unary systems. Each unary potential has 14 unary model parameters: cohesive energy (E_c), equilibrium nearest-neighbor distance (r_e), bulk modulus (B), a parameter (d) related to the pressure derivative of bulk modulus (∂B/∂P), decay lengths (β⁽⁰⁾, β⁽¹⁾, β⁽²⁾, β⁽³⁾), weight factors (t⁽¹⁾, t⁽²⁾, t⁽³⁾), a parameter (A) for embedding function, and many-body screening parameters (C_min, C_max). The binary model parameters are E_c, r_e, B, d, the ratio between atomic electron density scaling factors (ρ₀), four C_min, and four C_max. Detailed descriptions of the model parameters and the 2NN MEAM formalism are provided in Appendix A.

In this study, the 2NN MEAM potential for the Mg–Mn binary system was developed using the 2NN MEAM potentials for pure Mg¹⁸⁾ and Mn (Table 1),¹⁹⁾ which reproduce the fundamental material properties of their associated systems. The binary model parameters were determined to reproduce the fundamental material properties of the Mg–Mn system. Among the parameters, the E_c, r_e, B, and d parameter values constitute the universal equation of state for a selected reference structure. Thus, they can be determined from the corresponding experimental data for the reference structure if the structure is a stable compound in a phase diagram. The reference structure in the (2NN) MEAM formalism only can be a structure whose second nearest neighboring atoms are composed of the same element such as B1, B2, B3, and L1₂ structures. However, the Mg–Mn binary phase diagram does not contain any stable intermetallic compounds.⁶⁾ Thus, we searched for information on metastable Mg_xMn_y compounds with these structures and found a first-principles calculation that reports the cohesive energy, lattice parameters, and bulk modulus of a hypothetical L1₂ Mg₃Mn compound.²⁰⁾ Since the information was sufficient to determine the values of the E_c, r_e, and B parameters, we chose the L1₂ Mg₃Mn compound as the reference structure of the Mg–Mn potential.

Table 1 2NN MEAM potential parameter sets for pure Mg¹⁸⁾ and Mn.¹⁹⁾ The units of the cohesive energy E_c, equilibrium nearest-neighbor distance r_e, and bulk modulus B are eV, Å, and 10¹² dyne/cm², respectively. The reference structures are hcp for Mg and bcc for Mn.

The model parameters for the Mg–Mn system were determined using the reference structures. The r_e parameter value was optimized to fit the lattice parameter of the Mg-rich hcp solid solutions. The E_c and B parameters were assigned values obtained from the first-principles calculation for the cohesive energy and bulk modulus of the hypothetical L1₂ Mg₃Mn compound,²⁰⁾ respectively. The C_min and C_max parameter values were set from each assumption presented in Table 2. The value of the d parameter was set as the mean value of its constituent pure elements. The ρ₀^Mg:ρ₀^Mn ratio was given a default assumed value. The final potential parameter set determined for the Mg–Mn system is listed in Table 2.

Table 2 A 2NN MEAM potential parameter set for the Mg–Mn binary system. The units of the cohesive energy E_c, equilibrium nearest-neighbor distance r_e, and bulk modulus B are eV, Å, and 10¹² dyne/cm², respectively.

3. Calculations of Fundamental Material Properties

In this section, to examine the reliability and transferability of the presently developed potential, the fundamental material properties of the Mg–Mn binary system are calculated and compared with data from the literature. The reliability was examined through properties related to the fitting quality. As mentioned before, there were no stable intermetallic compounds in the Mg–Mn phase diagram,⁶⁾ so the fitting quality of the Mg–Mn system was examined by calculating the properties of stable solutions and hypothetical compounds. The transferability of the potential was examined to evaluate whether it could describe the dislocation behavior of the Mg–Mn alloy and could be used at finite temperatures.

Except for the enthalpy of mixing of liquid solutions and internal energy of compounds and solid solutions, the material properties were calculated at 0 K using in-house molecular statics and dynamics code, KISSMD,²¹⁾ and a large-scale atomic/molecular massively parallel simulator (LAMMPS) package.²²⁾ A periodic boundary condition was set in all directions, except for the generalized stacking fault energy (GSFE) calculation. The radial cutoff distance was set to 6.0 Å which is a larger value than the second nearest-neighbor distances of Mg and Mn.

3.1 Fitting quality of the Mg–Mn potential

To describe a system, an interatomic potential should reproduce the structural, elastic, and thermodynamic properties, which are the most fundamental material properties, of the stable phases in the system. In this respect, we first confirmed the reproducibility of the structural property of the presently developed Mg–Mn potential by calculating the lattice parameters of the Mg-rich hcp solid solutions owing to the lack of stable intermetallic compounds in the Mg–Mn system.⁶⁾ For this calculation, we used eight random, disordered hcp solid solution samples containing 0 to 2.8 at% Mn, and calculated their lattice parameters at 0 K. Each sample contained 4000 atoms. The calculation results are presented in Fig. 1. The calculated lattice parameters decrease in value with the Mn content on both axes. The present calculations on the a-axis were consistent with the experimental data,^23–25) whereas those on the c-axis were a little lower than the experimental data. Nevertheless, the decrement ratios of the calculated lattice parameters on both axes were in good agreement with the experimental data.

Fig. 1

Lattice parameters of hcp Mg–Mn alloys at 0 K according to the present 2NN MEAM potential, in comparison with experimental data.^23–25)

The reproducibility of the elastic property was confirmed by calculating the elastic moduli values of a hypothetical L1₂ Mg₃Mn intermetallic compound, and the values are listed in Table 3, in comparison with a first-principles calculation.²⁰⁾ The calculated bulk modulus was consistent with the first-principles calculation.

Table 3 Elastic moduli (GPa) of a hypothetical L1₂ Mg₃Mn intermetallic compound according to the present 2NN MEAM potential, in comparison with a first-principles calculation (F.P. calc.).

The reproducibility of the thermodynamic property was confirmed by testing whether the present Mg–Mn potential could reproduce the Mg–Mn system that lacked any stable intermetallic compounds but contained a liquid miscibility gap at high temperatures.⁶⁾ The absence of stable compounds was identified by calculating the enthalpies of formation of hcp and cbcc solid solutions and many hypothetical intermetallic compounds, and the values are presented in Fig. 2(a), in comparison with a first-principles calculation.²⁰⁾ The calculated enthalpy of formation of the hypothetical L1₂ Mg₃Mn intermetallic compound was consistent with the first-principles calculation.²⁰⁾ Moreover, the calculated enthalpies of formation of all the hypothetical compounds were positive, which meant that they were unstable under the present Mg–Mn potential. Although the hypothetical compounds could not contain all the compounds in the whole materials, they contained the representative intermetallic compounds frequently formed in Mg alloys. Therefore, the present calculation results might reproduce the Mg–Mn system without stable compounds. Similar to the hypothetical compounds, the calculated enthalpies of formation of the hcp and cbcc solid solutions were all positive. This meant that the matrix in the Mg–Mn alloys was not composed of one kind of solid solution, i.e., the hcp solid solution or the cbcc solid solution, which is in agreement with the phase diagram of the Mg–Mn system whose matrix consists of the Mg-rich hcp solid solution and the Mn-rich cbcc solid solution.⁶⁾ The enthalpies of mixing of liquid solutions were calculated to check the presence of the liquid miscibility gap using 13 random, disordered hcp solid solution samples containing 0–100 at% Mn. Each sample contained 4000 atoms. These solid samples were converted into liquid samples by raising their temperature to 2000 K. The liquid samples were cooled to 1600 K (as per temperatures in comparable literature).^6,26) The calculation results are presented in Fig. 2(b), in comparison with CALPHAD calculations.^6,26) While the general enthalpy of mixing shows parabolic curves, the calculated enthalpy of mixing appears to contain a straight line in the middle of the parabolic curve, which appears from 20 at% Mn to 90–95 at% Mn. This straight line is representative of a miscibility gab and also appears in the CALPHAD calculations.^6,26) To check whether a phase separation occurs in the Mg–Mn liquid phase, we observed the atomic configuration of the Mg–50 at% Mn alloy simulated from the Parrinello/Rahman NPT ensemble at 1600 K for 250 ps. Figure 3 depicts that post-simulation, Mg and Mn atoms were divided into the Mg-rich liquid solution and the Mn-rich liquid solution. Although there is a discrepancy between the calculated enthalpy of mixing and the CALPHAD calculations (Fig. 2(b)), the present Mg–Mn potential could reproduce the phase separation phenomenon and reasonably predict the miscibility gap section.

Fig. 2

(a) Enthalpies of formation of metastable hcp and cbcc solid solutions and hypothetical compounds at 0 K and (b) enthalpies of mixing of liquid solutions at finite temperatures of the Mg–Mn binary system according to the present 2NN MEAM potential, in comparison with a first-principles calculation²⁰⁾ and CALPHAD calculations.^6,26)

Fig. 3

Snapshots of (a) the initial and (b) final structures of the Mg–50 at% Mn alloy sample simulated from the Parrinello/Rahman NPT ensemble at 1600 K for 250 ps. Colors represent the element type of atoms (Green: Mg, Purple: Mn). The number of atoms in the sample is 4000.

3.2 Transferability of the developed Mg–Mn potential

In the previous section, we confirmed that the developed Mg–Mn potential can reproduce the structural, elastic, and thermodynamic properties of the stable phases in the Mg–Mn system. These properties are important for describing the alloy system but do not guarantee that the developed potential can reproduce the dislocation behavior of the Mg–Mn alloy. To confirm the reproducibility, we calculated the GSFE that is related to the dislocation behaviors of metallic alloys.²⁷⁾ The GSFE of pure Mg was also calculated for comparison. Among the stacking fault types in hcp metals, we investigated the I₂-type stacking fault (…ABABCACA…), which is directly created by shear deformation. The calculation was conducted for representative slip planes and directions in hcp metals: basal $\{ 0001\} \langle 10\bar{1}0\rangle $, basal $\{ 0001\} \langle 11\bar{2}0\rangle $, prismatic $\{ 10\bar{1}0\} \langle 11\bar{2}0\rangle $, pyramidal I $\{ 10\bar{1}1\} \langle 11\bar{2}0\rangle $, pyramidal I $\{ 10\bar{1}1\} \langle 11\bar{2}3\rangle $, and pyramidal II $\{ 11\bar{2}2\} \langle 11\bar{2}3\rangle $ slip systems. The basal, prismatic, pyramidal I, and pyramidal II samples consist of 48 layers of the (0001) plane, 40 layers of the $(10\bar{1}0)$ plane, 36 layers of the $(10\bar{1}1)$ plane, and 40 layers of the $(11\bar{2}2)$ plane, respectively. The samples for the alloy system were created by replacing several Mg atoms on the stacking fault plane with Mn atoms. The solute concentration on the SF plane was 25.0 at%, which is equivalent to that from first-principles calculation²⁸⁾ with many comparable data. The detailed calculation procedure is described in our previous works.^16,17) The calculated GSFE curves for the slip systems are shown in Fig. 4. The reliability of the calculated GSFE was verified by comparing their stable and unstable stacking fault energies (SFEs) with first-principles calculations^28–32) (Table 4). The calculated stable SFE of the basal $\{ 0001\} \langle 10\bar{1}0\rangle $ slip system is larger in the order of pure Mg and Mg–Mn. The SFE order is consistent with those from first-principles calculations^28,29) but is opposite to that of a first-principles calculation.³⁰⁾ The calculated stable SFE of the pyramidal II $\{ 11\bar{2}2\} \langle 11\bar{2}3\rangle $ slip system is larger in the order of pure Mg and Mg–Mn. Though there is no comparable data to confirm the consistency of this SFE order, these values are similar to those from a first-principles calculation.³¹⁾ The calculated unstable SFEs of the basal $\{ 0001\} \langle 10\bar{1}0\rangle $ and $\{ 0001\} \langle 11\bar{2}0\rangle $ slip systems are both smaller in the order of pure Mg and Mg–Mn, which is opposite to the orders from the first-principles calculations.^28,29) The calculated unstable SFE of the prismatic $\{ 10\bar{1}0\} \langle 11\bar{2}0\rangle $ slip system is larger in the order of pure Mg and Mg–Mn, which coincides with the SFE order from a first-principles calculation.²⁸⁾ The calculated unstable SFEs of the pyramidal I $\{ 10\bar{1}1\} \langle 11\bar{2}0\rangle $ and pyramidal I $\{ 10\bar{1}1\} \langle 11\bar{2}3\rangle $ slip systems are both larger in the order of pure Mg and Mg–Mn, and those of the pyramidal II $\{ 11\bar{2}2\} \langle 11\bar{2}3\rangle $ slip system are smaller in the order of pure Mg and Mg–Mn. The consistency of those SFE orders is difficult to confirm owing to the lack of comparable data, but those SEF values are similar to those from first-principle calculations.^31,32)

Fig. 4

Calculated GSFE (mJ/m²) curves of pure Mg and Mg–Mn, models on the (a) basal plane along the $\langle 10\bar{1}0\rangle $ direction, (b) basal plane along the $\langle 11\bar{2}0\rangle $ direction, (c) prismatic plane along the $\langle 11\bar{2}0\rangle $ direction, (d) pyramidal I plane along the $\langle 11\bar{2}0\rangle $ direction, (e) pyramidal I plane along the $\langle 11\bar{2}3\rangle $ direction, and (f) pyramidal II plane along the $\langle 11\bar{2}3\rangle $ direction, using the present 2NN MEAM potential.

Table 4 Stable and unstable SFE (mJ/m²) of pure Mg and Mg–Mn models for the basal, prismatic, pyramidal I, and pyramidal II slip systems according to the present 2NN MEAM potentials, in comparison with first-principles calculations.

In addition to the GSFE, we confirmed whether the presently developed Mg–Mn potential is reliable at finite temperatures for application in various future simulations. The reliability for the liquid phases was identified in the previous section (Fig. 2(b)), but that for the solid phases was not confirmed. To identify the reliability for the solid phases, we investigated whether the Mg_xMn_y hypothetical compounds and metastable solid solutions were unstable at finite temperatures. The stability of the compounds and solid solutions was identified using their internal energy after heating and rapid cooling to 0 K. Notably, the internal energy of the Mg_xMn_y intermetallic compounds and solid solutions should be less stable than the convex hull of pure Mg and pure Mn at the given temperatures and compositions. To confirm this, we heated the Mg_xMn_y hypothetical compound samples and the metastable hcp and cbcc solid solution samples to 1000 K and rapidly cooled them to 0 K from each heating temperature. The samples contained approximately 2000 atoms and the heating interval was 100 K. At each heating temperature, each sample was equilibrated for 30 ps. As shown in Fig. 5 and Fig. 6, the calculated internal energy of the whole Mg_xMn_y hypothetical compounds and metastable solid solutions was less negative than the convex hull, meaning that they were all unstable from 0 to 1000 K.

Fig. 5

Change in the internal energy of stable (filled symbols) and metastable (open symbols) phases of (a) Mg₃Mn, (b) Mg₂Mn, (c) MgMn, (d) MgMn₂, and (e) MgMn₃ compounds after heating and rapid cooling to 0 K from each heating temperature.

Fig. 6

Change in the internal energy of metastable hcp and cbcc solid solutions and hypothetical compounds in (a) Mg₉Mn, (b) Mg₈Mn₂, (c) Mg₇Mn₃, (d) Mg₆Mn₄, (e) Mg₅Mn₅, (f) Mg₄Mn₆, (g) Mg₃Mn₇, (h) Mg₂Mn₈, and (i) MgMn₉ compositions after heating and rapid cooling to 0 K from each temperature.

Taken together, the presently developed Mg–Mn interatomic potential can reproduce the structural, elastic, and thermodynamic properties of (meta)stable phases in the Mg–Mn system, and reasonably reproduce the GSFE of the Mg–Mn alloy. Furthermore, this potential can be used not only at 0 K but also at finite temperatures up to 1000 K. Therefore, the presently developed Mg–Mn interatomic potential can be utilized for studying the dislocation behavior of the Mg–Mn alloy, which will contribute to revealing the effect of Mn on the RT formability of Mg.

4. Conclusion

Based on the 2NN MEAM formalism, an interatomic potential for the Mg–Mn binary system was developed to infer the effect of Mn on the RT formability of Mg. The developed potential reasonably reproduced the structural, elastic, and thermodynamic properties of the compound and solution phases of the associated alloy system. This potential also described the GSFE of the basal, prismatic, pyramidal I, and pyramidal II slip systems of the Mg–Mn alloy, in reasonable agreement with first-principles calculations. Furthermore, it can be utilized from 0 K to 1000 K. The developed potential could be utilized to investigate the dislocation behavior of the Mg–Mn alloy and to elucidate the effect of Mn on the formability of Mg, which will further contribute to the enhancement of the RT formability and strength of Mg alloys.

Acknowledgments

This work was supported by the National Research Foundation of Korea (NRF) grants funded by the Korea government (MSIP) (No. 2021M3A7C2089767, 2021M3H4A6A01045764) and the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Science, ICT & Future Planning (No. 2016R1A2B4006680).

REFERENCES

• 1)   J.J.  Kim and  D.S.  Han: Mater. Trans. 49 (2008) 894–897. doi:10.2320/matertrans.MC200731
• 2)   J.  Hirsch: Mater. Trans. 52 (2011) 818–824. doi:10.2320/matertrans.L-MZ201132
• 3)   H.  Haferkamp,  R.  Boehm,  U.  Holzkamp,  C.  Jaschik,  V.  Kaese and  M.  Niemeyer: Mater. Trans. 42 (2001) 1160–1166. doi:10.2320/matertrans.42.1160
• 4)  A.A. Luo, M.O. Pekguleryuz, K.U. Kainer and A.A. Kaya: Fundam. Magnes. Alloy Metall., (Woodhead Publishing, Cambridge, 2013) pp. 266–316.
• 5)  M.O. Pekguleryuz, K.U. Kainer and A.A. Kaya: Fundamentals of Magnesium Alloy Metallurgy, (Woodhead Publishing, Cambridge, 2013).
• 6)   J.  Gröbner,  D.  Mirkovic,  M.  Ohno and  R.  Schmid-Fetzer: J. Phase Equilibria Diffus. 26 (2005) 234–239. doi:10.1007/s11669-005-0110-8
• 7)   A.  Shukla and  A.D.  Pelton: J. Phase Equilibria Diffus. 30 (2009) 28–39. doi:10.1007/s11669-008-9426-5
• 8)   P.  Ghosh and  M.  Medraj: Calphad 41 (2013) 89–107. doi:10.1016/j.calphad.2013.01.008
• 9)   H.  Somekawa,  A.  Kinoshita and  A.  Kato: Mater. Sci. Eng. A 732 (2018) 21–28. doi:10.1016/j.msea.2018.06.098
• 10)   X.  Gu,  Y.  Zheng,  Y.  Cheng,  S.  Zhong and  T.  Xi: Biomaterials 30 (2009) 484–498. doi:10.1016/j.biomaterials.2008.10.021
• 11)   J.  Wang,  S.K.  Yadav,  J.P.  Hirth,  C.N.  Tomé and  I.J.  Beyerlein: Mater. Res. Lett. 1 (2013) 126–132. doi:10.1080/21663831.2013.792019
• 12)   H.  Feng,  Q.H.  Fang,  B.  Liu,  Y.  Liu,  Y.W.  Liu and  P.H.  Wen: Mech. Mater. 109 (2017) 26–33. doi:10.1016/j.mechmat.2017.03.015
• 13)  F.J. Humphreys and M. Hatherly: Recrystallization and Related Annealing Phenomena, Second ed., (Elsevier Ltd, Oxford, 2004).
• 14)   B.-J.  Lee and  M.I.  Baskes: Phys. Rev. B 62 (2000) 8564–8567. doi:10.1103/PhysRevB.62.8564
• 15)   H.-S.  Jang,  J.-K.  Lee,  A.J.S.F.  Tapia,  N.J.  Kim and  B.-J.  Lee: J. Magnes. Alloy. 10 (2022) 585–597. doi:10.1016/j.jma.2021.03.007
• 16)   H.-S.  Jang,  D.  Seol and  B.-J.  Lee: J. Alloy. Compd. 894 (2022) 162539. doi:10.1016/j.jallcom.2021.162539
• 17)   H.-S.  Jang,  D.  Seol and  B.-J.  Lee: J. Magnes. Alloy. 9 (2021) 317–335. doi:10.1016/j.jma.2020.09.006
• 18)   Y.M.  Kim,  N.J.  Kim and  B.-J.  Lee: Calphad 33 (2009) 650–657. doi:10.1016/j.calphad.2009.07.004
• 19)   Y.M.  Kim,  Y.H.  Shin and  B.J.  Lee: Acta Mater. 57 (2009) 474–482. doi:10.1016/j.actamat.2008.09.031
• 20)  https://materialsproject.org/materials/mp-978286/ (accessed Nov. 17, 2017).
• 21)  https://cmse.postech.ac.kr/home_2nnmeam (accessed March 20, 2020).
• 22)   S.  Plimpton: J. Comput. Phys. 117 (1995) 1–19. doi:10.1006/jcph.1995.1039
• 23)   E.  Schmid and  G.  Siebel: Metallwirtschaft 10 (1931) 923–925.
• 24)   R.S.  Busk: JOM 2 (1950) 1460–1464. doi:10.1007/BF03399173
• 25)   N.I.  Varich and  B.N.  Litvin: Phys. Met. Met. 16 (1963) 29–32.
• 26)  J. Tibballs, I. Ansara, A. Dinsdale and M. Rand: COST507-Thermochemical Database Light Met. Alloy., (European Commission EUR 18499, 1998) pp. 215–217.
• 27)   V.  Vítek: Philos. Mag. 18 (1968) 773–786. doi:10.1080/14786436808227500
• 28)   M.  Muzyk,  Z.  Pakiela and  K.J.  Kurzydlowski: Scr. Mater. 66 (2012) 219–222. doi:10.1016/j.scriptamat.2011.10.038
• 29)   S.L.  Shang,  W.Y.  Wang,  B.C.  Zhou,  Y.  Wang,  K.A.  Darling,  L.J.  Kecskes,  S.N.  Mathaudhu and  Z.K.  Liu: Acta Mater. 67 (2014) 168–180. doi:10.1016/j.actamat.2013.12.019
• 30)   W.Y.  Wang,  S.L.  Shang,  Y.  Wang,  Z.G.  Mei,  K.A.  Darling,  L.J.  Kecskes,  S.N.  Mathaudhu,  X.D.  Hui and  Z.K.  Liu: Mater. Res. Lett. 2 (2014) 29–36. doi:10.1080/21663831.2013.858085
• 31)   C.  Wang,  H.-Y.  Wang,  H.-Y.  Zhang,  X.-L.  Nan,  E.-S.  Xue and  Q.-C.  Jiang: J. Alloy. Compd. 575 (2013) 423–433. doi:10.1016/j.jallcom.2013.05.187
• 32)   T.  Hase,  T.  Ohtagaki,  M.  Yamaguchi,  N.  Ikeo and  T.  Mukai: Acta Mater. 104 (2016) 283–294. doi:10.1016/j.actamat.2015.11.045
• 33)   J.H.  Rose,  J.R.  Smith,  F.  Guinea and  J.  Ferrante: Phys. Rev. B 29 (1984) 2963–2969. doi:10.1103/PhysRevB.29.2963
• 34)   M.I.  Baskes: Phys. Rev. B 46 (1992) 2727–2742. doi:10.1103/PhysRevB.46.2727
• 35)   M.I.  Baskes: Mater. Chem. Phys. 50 (1997) 152–158. doi:10.1016/S0254-0584(97)80252-0
• 36)   B.-J.  Lee,  M.I.  Baskes,  H.  Kim and  Y.K.  Cho: Phys. Rev. B 64 (2001) 184102. doi:10.1103/PhysRevB.64.184102
• 37)   B.-J.  Lee,  J.-H.  Shim and  M.I.  Baskes: Phys. Rev. B 68 (2003) 144112. doi:10.1103/PhysRevB.68.144112
• 38)   B.-J.  Lee,  W.-S.  Ko,  H.-K.  Kim and  E.-H.  Kim: Calphad 34 (2010) 510–522. doi:10.1016/j.calphad.2010.10.007
• 39)   M.I.  Baskes: Mater. Sci. Eng. A 261 (1999) 165–168. doi:10.1016/S0921-5093(98)01062-4

Appendix  A

In the MEAM, the total energy of a system is approximated as   

\begin{equation} E = \sum_{i} \left[F_{i}(\bar{\rho}_{i}) + \frac{1}{2}\sum_{j(\neq i)} S_{ij}\phi_{ij}(R_{ij})\right] \end{equation} (A.1)

where F_i is the embedding function, $\bar{\rho }_{i}$ is the background electron density at site i. S_ij and ϕ_ij(R_ij) are the screening factor and the pair interaction between atoms i and j separated by a distance R_ij, respectively. For general calculations of energy, the functional forms for the two terms on the right-hand side of eq. (A.1), F_i and ϕ_ij, should be given. The background electron density at a site is computed considering directionality in bonding. A specific form is given to the embedding function F_i, but not to the pair interaction ϕ_ij. Instead a reference structure where individual atoms are on the exact lattice points is defined and the total energy per atom of the reference structure is estimated from the zero-temperature universal equation of state of Rose et al.³³⁾ Then, the value of the pair interaction is evaluated from the known values of the total energy per atom and the embedding energy, as a function of nearest-neighbor distance. In the original MEAM,³⁴⁾ only first nearest-neighbor interactions are considered. The neglect of the second nearest-neighbor and more distant nearest-neighbor interactions is made effective by the use of a strong many-body screening function.³⁵⁾ The consideration of the second nearest-neighbor interactions in the modified formalism (the second nearest-neighbor MEAM^14,36–38) is affected by adjusting the screening parameters so that the many-body screening becomes less severe.

It has been shown in eq. (A.1) that the energy in the MEAM is composed of two terms, the embedding function term and the pair interaction term. The embedding function is given the following form,   

\begin{equation} F(\bar{\rho}) = AE_{c}\frac{\bar{\rho}}{\bar{\rho}^{o}}\ln \frac{\bar{\rho}}{\bar{\rho}^{o}} \end{equation} (A.2)

where A is an adjustable parameter, E_c is the cohesive energy and $\bar{\rho }^{o}$ is the background electron density for a reference structure. The reference structure is a structure where individual atoms are on the exact lattice points. Normally, the equilibrium structure is taken as the reference structure for elements. The background electron density $\bar{\rho }_{i}$ is composed of spherically symmetric partial electron density, ρ_i⁽⁰⁾, and angular contributions, ρ_i⁽¹⁾, ρ_i⁽²⁾, ρ_i⁽³⁾. Each partial electron density term has the following form,^34,39)   

\begin{equation} (\rho_{i}^{(0)})^{2} = \left[ \sum_{j\neq i} S_{ij}\rho_{j}^{a(0)}(R_{ij})\right]^{2} \end{equation} (A.3a)

  

\begin{equation} (\rho_{i}^{(1)})^{2} = \sum_{\alpha} \left[\sum_{j\neq i} \frac{R_{ij}^{\alpha}}{R_{ij}}S_{ij}\rho_{j}^{a(1)}(R_{ij})\right]^{2} \end{equation} (A.3b)

  

\begin{align} (\rho_{i}^{(2)})^{2} &= \sum_{\alpha,\beta} \left[ \sum_{j \neq i} \frac{R_{ij}^{\alpha} R_{ij}^{\beta}}{R_{ij}^{2}}S_{ij}\rho_{j}^{a(2)}(R_{ij}) \right]^{2}{}\\ &\quad - \frac{1}{3}\left[\sum_{j\neq i} S_{ij}\rho_{j}^{a(2)}(R_{ij})\right]^{2} \end{align} (A.3c)

  

\begin{align} (\rho_{i}^{(3)})^{2}& = \sum_{\alpha,\beta,\gamma} \left[\sum_{j\neq i} \frac{R_{ij}^{\alpha} R_{ij}^{\beta} R_{ij}^{\gamma}}{R_{ij}^{3}}S_{ij}\rho_{j}^{a(3)}(R_{ij})\right]^{2} \\ &\quad - \frac{3}{5}\sum_{\alpha} \left[\sum_{j\neq i} \frac{R_{ij}^{\alpha}}{R_{ij}}{S_{ij}}\rho_{j}^{a(3)}(R_{ij})\right]^{2} \end{align} (A.3d)

Here, ρ_j^a(h) represent atomic electron densities from j atom at a distance R_ij from site i. R_ij^α is the α component of the distance vector between atoms j and i (α = x, y, z). The expression for (ρ_i⁽³⁾)² (eq. (A.3d)) is that modified later in order to make the partial electron densities orthogonal (eq. (A.4)). The way of combining the partial electron densities to give the total background electron density is not unique, and several expressions have been proposed.³⁵⁾ Among them, the following form that can be widely used without numerical error is taken in the present work.   

\begin{equation} \bar{\rho}_{i} = \rho_{i}^{(0)}G(\Gamma) \end{equation} (A.4)

where   

\begin{equation} G(\Gamma) = \frac{2}{1 + {e^{-\Gamma}}} \end{equation} (A.5)

and   

\begin{equation} \Gamma = \sum_{h=1}^{3} t_{i}^{(h)} \left[\frac{\rho_{i}^{(h)}}{\rho_{i}^{(0)}} \right]^{2} \end{equation} (A.6)

t_i^(h) are adjustable parameters. The atomic electron density is given as,   

\begin{equation} \rho_{j}^{a(h)}(R) = \rho_{0}e^{-\beta^{(h)}(R/r_{e} - 1)} \end{equation} (A.7)

where ρ₀ is a scaling factor, β^(h) are adjustable parameters and r_e is the nearest-neighbor distance in the equilibrium reference structure. The scaling factor doesn’t give any effect on calculations for elements, but can give effect on calculations for alloy systems.

In the MEAM no specific functional expression is given directly to ϕ(R). Instead, the atomic energy (total energy per atom) is evaluated by some means as a function of nearest-neighbor distance. Then, the value of ϕ(R) is computed from known values of the total energy and the embedding energy, as a function of nearest-neighbor distance.

Let’s consider a reference structure once again. Here, every atom has the same environment and the same energy. If up to second nearest-neighbor interactions are considered as is done in the 2NN MEAM,^14,36) the total energy per atom in a reference structure can be written as follows:   

\begin{equation} E^{a}(R) = F(\bar{\rho}^{o}(R)) + \frac{Z_{1}}{2}\phi (R) + \frac{Z_{2}S}{2}\phi (aR) \end{equation} (A.8)

where Z₁ and Z₂ are the number of first and second nearest-neighbor atoms, respectively. S is the screening factor for second nearest-neighbor interactions (the screening factor for first nearest-neighbor interactions is 1), and a is the ratio between the second and first nearest-neighbor distances. It should be noted that for a given reference structure S and a are constants, and the total energy and the embedding energy become functions of only nearest-neighbor distance R. On the other hand, the energy per atom for a reference structure can be obtained from the zero-temperature universal equation of state of Rose et al.³³⁾ as a function of nearest-neighbor distance R.   

\begin{equation} E^{u}(R) = -E_{c}(1 + a^{*} + da^{*3})e^{-a^{*}} \end{equation} (A.9)

where d is an adjustable parameter, and   

\begin{equation} a^{*} = \alpha (R/r_{e} - 1) \end{equation} (A.10)

and   

\begin{equation} \alpha = \left(\frac{9B\Omega}{E_{c}}\right)^{1/2}. \end{equation} (A.11)

E^u(R) is the universal function for a uniform expansion or contraction in the reference structure, B is the bulk modulus and Ω is the equilibrium atomic volume.

Basically, the pair potential between two atoms separated by a distance R, ϕ(R), can be obtained by equating eq. (A.8) and eq. (A.9). However, it is not trivial because eq. (A.8) contains two pair potential terms. In order to derive an expression for the pair interaction, ϕ(R), another pair potential, ψ(R), is introduced:   

\begin{equation} E^{u}(R) = F(\bar{\rho}^{o}(R)) + \frac{Z_{1}}{2}\psi (R) \end{equation} (A.12)

where   

\begin{equation} \psi (R) = \phi (R) + \frac{Z_{2}S}{Z_{1}}\phi (aR). \end{equation} (A.13)

ψ(R) can be computed from eq. (A.11) as a function of R, as follows:   

\begin{equation} \psi (R) = \frac{2}{Z_{1}}[E^{u}(R) - F(\bar{\rho}^{o}(R))], \end{equation} (A.14)

and the expression for the pair potential ϕ(R) is obtained from eq. (A.12) as follows:   

\begin{equation} \phi (R) = \psi (R) + \sum_{n=1} (-1)^{n} \left(\frac{Z_{2}S}{Z_{1}} \right)^{n}\psi (a^{n}R). \end{equation} (A.15)

Here, the summation is performed until a correct value of atomic energy is obtained for the equilibrium reference structure.

It should be noted here that the original first nearest neighbor MEAM is a special case (S = 0) of the present 2NN MEAM. In the original MEAM, the neglect of the second nearest-neighbor interactions is made by the use of a strong many-body screening function.³⁵⁾ In the same way, the consideration of the second nearest-neighbor interactions in the modified formalism is affected by adjusting the many-body screening function so that it becomes less severe. In the MEAM, the many-body screening function between atoms i and j, S_ij, is defined as the product of the screening factors, S_ikj, due to all other neighbor atoms k:   

\begin{equation} S_{ij} = \prod_{k\neq i,j} {S_{ikj}} \end{equation} (A.16)

The screening factor S_ikj is computed using a simple geometric construction. Imagine an ellipse on an x, y plane, passing through atoms, i, k and j with the x-axis of the ellipse determined by atoms i and j. The equation of the ellipse is given by,   

\begin{equation} x^{2} + \frac{1}{C}y^{2} = \left(\frac{1}{2}R_{ij} \right)^{2}. \end{equation} (A.17)

For each k atom, the value of parameter C can be computed from relative distances among the three atoms, i, j and k, as follows:   

\begin{equation} C = \frac{2(X_{ik} + X_{kj}) - (X_{ik} - X_{kj})^{2} - 1}{1 - (X_{ik} - X_{kj})^{2}} \end{equation} (A.18)

where X_ik = (R_ik/R_ij)² and X_kj = (R_kj/R_ij)². The screening factor, S_ikj is defined as a function of C as follows:   

\begin{equation} S_{ikj} = f_{c}\left[\frac{C - C_{\min}}{C_{\max} - C_{\min}}\right] \end{equation} (A.19)

where C_min and C_max are the limiting values of C determining the extent of screening and the smooth cutoff function is   

\begin{equation} \begin{array}{ll} f_{c}(x) = 1 & x \geq 1\\ {[}1 - (1-x)^{4}]^{2} & 0 < x < 1\\ 0 & x \leq 0. \end{array} \end{equation} (A.20)

The basic idea for the screening is as follows: First define two limiting values, C_max and C_min (C_max > C_min). Then, if the atom k is outside of the ellipse defined by C_max, it is thought that the atom k does not have any effect on the interaction between atoms i and j. If the atom k is inside of the ellipse defined by C_min it is thought that the atom k completely screens the i-j interaction, and between C_max and C_min the screening changes gradually. In the numerical procedure of simulation the electron density and pair potential are multiplied by the screening function S_ij, as is done in eq. (A.1) and eqs. (A.3a)–(A.3d). Therefore, S_ij = 1 and S_ij = 0 mean that the interaction between atoms i and j is unscreened and completely screened, respectively. In addition to the many-body screening function, a radial cutoff function which is given by f_c[(r_c − r)/Δr] where r_c is the cutoff distance and Δr (0.1 Å) is the cutoff region, is also applied to the atomic electron density and pair potential.³⁵⁾ The radial cutoff distance is chosen so that it doesn’t have any effect on the calculation results due to the many-body screening. This is only for the computational convenience, that is, to save the computation time.

To describe a binary alloy system, in addition to the descriptions for individual elements, the pair interaction between different elements should be determined. For this, a similar technique used to determine the pair interaction for pure elements is applied. A perfectly ordered binary intermetallic compound, where one type of atom has only same type of atoms as second nearest-neighbors, is considered as a reference structure. The Bl (NaCl type) or B2 (CsCl type) ordered structure can be a good example. For such a reference structure, the total energy per atoms (for l/2 i atom + 1/2 j atom) is given as follows:   

\begin{align} E_{ij}^{u}(R) &= \frac{1}{2}\bigg[ F_{i}(\bar{\rho}_{i}) + F_{j}(\bar{\rho}_{j}) + Z_{1}\phi_{ij}(R) \\ &\quad + \frac{Z_{2}}{2}( S_{i}{\phi_{ii}}(aR) + S_{j}\phi_{jj}(aR)) \bigg] \end{align} (A.21)

where Z₁ and Z₂ are the numbers of first and second nearest-neighbors in the reference structure, respectively. S_i and S_j are the screening functions for the second nearest-neighbor interactions between i atoms and between j atoms, respectively, and a is the ratio between the second and first nearest-neighbor distances in the reference structure. The pair interaction between i and j can now be obtained in the following form:   

\begin{align} \phi_{ij}(R) &= \frac{1}{Z_{1}}\bigg[2E_{ij}^{u}(R) - F_{i}(\bar{\rho}_{i}) - F_{j}(\bar{\rho}_{j}) \\ &\quad - \frac{Z_{2}}{2}(S_{i}\phi_{ii}(aR) + S_{j}\phi_{jj}(aR))\bigg]. \end{align} (A.22)

The embedding functions F_i and F_j can be readily computed. The pair interactions ϕ_ii and ϕ_jj between the same types of atoms can also be computed from descriptions of individual elements. To obtain $E_{ij}^{u}(R)$, the universal equation of state, eqs. (A.9)–(A.11), is considered once again for the binary reference structure. In this case, E_c, r_e and B correspond to the cohesive energy, equilibrium nearest-neighbor distance and bulk modulus of the binary reference structure.

In addition to the parameters for the universal equation of state, two more model parameter groups, C_min and C_max, must be determined to describe alloy systems. Each element has its own value of C_min and C_max. C_min and C_max determine the extent of screening of an atom (k) to the interaction between two neighboring atoms (i and j). For pure elements, the three atoms are all the same type (i-j-k = A-A-A or B-B-B). However, in the case of binary alloys, one of the interacting atoms and/or the screening atoms can be different types (there are four cases: i-k-j = A-B-A, B-A-B, A-A-B and A-B-B). Different C_min and C_max values may have to be given in each case. Another, the last model parameter necessary for binary alloy systems is the atomic electron density scaling factor ρ_o (see eq. (A.7)). For an equilibrium reference structure (R = r_e), the values of all atomic electron densities become ρ_o. This is an arbitrary value and does not have any effect on calculations for pure elements. This parameter is often omitted when describing the potential model for pure elements. However, for alloy systems, especially for systems where the constituent elements have different coordination numbers, the scaling factor (the ratio of the two values) has a great effect on calculations.