Evaluating the applicability of classical and neural network interatomic potentials for modeling body centered cubic polymorph of magnesium

Magnesium (Mg) is one of the most abundant metallic elements in nature and presents attractive mechanical properties in the industry. Particularly, it has a low density and relatively high strength/weight and stiffness/weight ratios, which make it one of the most attractive lightweight metals. However, the huge potential of Mg is restricted by its low ductility, associated with its hexagonal close packed (hcp) structure. This problem can be solved if Mg adopts the body centered cubic (bcc) structure, which is stable at high pressure or in confinement with stiff bcc metals like Nb. Molecular dynamics method is a magnificent tool to study material’s structure and deformation mechanisms at the atomic level, however, requiring accurate interatomic potentials. The majority of the interatomic potentials available in the literature for Mg have only been fitted to the properties of its stable hcp phase. In the present work, we perform systematic study of applicability of currently available Mg potentials to modeling the properties of metastable bcc polymorph of Mg, taking into account cohesive energy curves, elastic constants, stacking fault energies, and phonon dispersion curves. We conclude that the modified embedded atom method (MEAM) potentials are the most suitable for investigating bcc Mg in Mg/Nb nano-composites, while the properties of high-pressure bcc Mg would be better modeled by neural network interatomic potentials after different local atomic environments corresponding to bcc Mg being included into the fitting database.


Introduction
Magnesium (Mg) is one of the most commonly used structural metals due to its abundance in nature and characteristic mechanical properties, especially its low density (1.74 g cm −3 ) [1]. Its excellent strength-to-weight and stiffness-to-weight ratios and large power storage capacity make Mg a promising candidate for aerospace and automotive applications, degradable implant materials and energy storage devices [2]. Mg is the lightest metal in nature: a 33% lighter than aluminum (Al), 60% lighter than titanium (Ti), and 75% lighter than steel. Consequently, replacing these metals in the aerospace and automotive industries can significantly reduce the weight of these vehicles, which results in reduced fuel consumption and CO 2 emission. However, Mg presents a relatively low melting temperature and a high brittleness with respect to other structural metals, so its range of application is limited. These properties can be slightly improved by alloying the material with other elements, including Al, Zn or Ti, but further research is still necessary to develop cheap, environment-friendly and efficient structural materials. The low ductility observed in Mg and Mg-based alloys is associated with its hexagonal close-packed (hcp) lattice structure, with few easily-activated slip systems [2,3], which inhibit the formation and slip of dislocations with a resulting strain hardening [4,5]. The active slip modes in Mg are largely confined to its close-packed basal planes [6,7], meaning that strains applied out of the basal plane cannot be supported by such dislocations [8]. Other slip systems (prismatic and pyramidal) exist but are considerably harder to activate, leading to the low ductility/formability of Mg. Deformation twinning is an important alternative deformation mode that contributes to the plasticity of Mg, and {1012} tension twins are profusely observed during plastic deformation of Mg [9,10]. Such tension twins can grow to larger sizes and consume the entire parent structure [11]. Ultimately, the intricate competition between slip and twinning is crucial for understanding and controlling the mechanical properties of Mg. However, further research is still necessary to improve the mechanical properties of hcp Mg.
In the presence of interface strains or under extreme external loads, the crystal structure of Mg can be transformed and stabilized from its natural hcp to bcc crystal structure [12,13]. This bcc phase exhibits significantly higher yield stresses and strains-to-failure, especially at ambient temperature and pressure when Mg forms a nanolaminate architecture with another bcc metal, such as Nb. Therefore, the use of Mg in its bcc phase is a promising alternative for the development of improved structural materials. Mg is the only alkaline earth metal that does not present the hcp → bcc phase transition with increasing temperature at ambient pressure and, in the pure material, this phase transition is only observed with increasing pressure [12]. At room temperature, the phase transition from hcp to bcc occurs at ≈50 GPa. The hcp-bcc phase line depends strongly on temperature and there is a triple point (bcc, hcp, liquid) at ≈4 GPa and ≈1200 K [14], while the melting temperature of Mg at zero pressure is 922 K [15]. Pure bcc Mg is a 50% stronger than pure hcp Mg and can resist up to a 60% larger strains to failure [16], which makes bcc Mg a highly attractive pseudomorphic phase. At ambient temperature and pressure, this bcc structure is stable in thin layers, such as in thin Mg/Nb multilayers. In this system, the bcc → hcp phase transition occurs if the layer thickness is larger than 6 nm [17]. If the thickness is between 7 and 10 nm, both phases, hcp and bcc, coexist, and Mg presents the hcp structure if the thickness is larger than 10 nm, with a thin bcc Mg layer (2-3 atomic planes) present at the Mg/Nb interface [17]. The stability of this bcc structure is due to interface strains and short-range interactions with bcc Nb. However, the origin and atomic-scale mechanisms behind this phase transition are still not well understood and the structural behaviour of the pure bcc Mg phase is still not known due to the difficulties in isolating and measuring its mechanical properties.
Computer simulations, including molecular dynamics (MD) simulations, density functional theory (DFT) calculations, or phase-field (PF) modeling, can be employed as a tool to elucidate the phenomena behind the bcc → hcp phase transition in Mg and the stabilization of the bcc-Mg phase at different length scales. The DFT method has been used to characterize the elastic properties of pure bcc Mg through the calculation of the elastic constants and the moduli of the perfect bcc crystal [16,18,19]. Although DFT calculations and ab initio MD simulations provide the best predictive accuracy of mechanical, thermal, electrical, and transport properties of materials [20], they are computationally demanding and, thus, limited to small simulation cells in order of a few nanometers and short simulation time in order of picoseconds. As a consequence, the study of the time evolution of the bcc-Mg phase in thin and thick Mg/Nb multilayers and in pure Mg at large pressure is not possible with the DFT method. On the other hand, PF simulations require the previous knowledge and understanding of all mechanisms at the atomic-and mesoscale, and do not provide additional information about the phenomena occurring at the atomic level. Finally, classical MD simulations have been proved to be an excellent tool for the characterization of mechanical, thermal, and transport properties of materials at the nanoscale, allowing for simulations with millions or even billions of atoms over the range of nano-or even microseconds on modern computing infrastructures [21][22][23][24][25]. Thus, MD simulations are an excellent tool to analyze the atomic-scale processes that cannot be understood directly from experimental observations. In MD simulations, all atoms are considered as interacting particles and Newton's equations of motion are integrated numerically to study time evolution of the system [26]. Mechanical, transport and thermal properties can be calculated directly from these trajectories using many advancements of statistical physics [26]. The interatomic interactions are described using interatomic potentials, i.e. mathematical functions that capture the behavior of the different atoms in the simulation box through the interaction of each atom with its neighbors. Many empirical functional forms have been proposed in the literature to model interatomic interactions and chemical bonds, with model parameters that are commonly fitted to reproduce experimental data or first-principles calculations [27,28]. The recent emergence of machine learning/neural network potentials promises the DFT accuracy at the fraction of computing time by being 3-5 orders of magnitude faster, however demonstrating much slower (1-2 orders of magnitude) performance than classical/empirical potentials [29]. Since the appearance of these machine learning potentials, the predictive power for complex long-range and time-averaged properties of materials has become feasible to evaluate, overcoming the limitations of ab initio and classical atomistic modeling methods.
MD simulations have been widely used in the study of different phenomena in Mg and Mg-based alloys, including twin nucleation [10], twin growth [30,31], fracture [32], interactions with dislocations [33,34], formation of grain boundaries [35], crystal plasticity [32,36] or alloying [37]. However, the majority of MD simulations with Mg have been performed to study the hcp phase, and the interatomic potentials available in the literature have not been fitted to properties of the bcc phase. In the present work, we analyze the capabilities of the Mg interatomic potentials from the NIST interatomic potentials repository [38], together with  2 Zhou et al [43] EAM 3 Sun et al [44] EAM 4 Wilson et al [45] EAM 5 Pei et al [46] MEAM 1 Kim et al [47] MEAM 2 Dickel et al [48] MEAM 3 Ahmad et al [49] MEAM 4 Ahmad et al [49] TB-SMA 1 Cleri et al [50] TB-SMA 2 Li et al [51] NNp 1 Stricker et al [39] NNp 2 Dickel et al [40] the emerging new machine learning/neural network potentials, in the prediction of thermodynamic, mechanical, and transport properties of bcc Mg, presenting their advantages and flaws. The interatomic potentials used in the present work are introduced in section 2 and tested in section 3.

Methods
In MD simulations, interatomic potentials define the interactions between atoms based on their relative positions in the system, thus giving the expression for the total potential energy in the system, from which all the forces acting on each atom at each moment of time can be derived.
In the present work, the structural properties of bcc Mg are studied with fourteen interatomic potentials available in the literature. This set of potentials includes potentials with different functional forms and different parameterizations for the same functional form: five embeddedatom method (EAM) potentials, four modified embedded-atom method (MEAM) potentials, one angular dependent potential (ADp), two neural network potentials (NNp) and two tightbinding, with the second-moment approximation, (TB-SMA) models. These potentials have been proved to reproduce the mechanical properties of hcp Mg, and, in the case of the NNps, they could also reproduce the fundamental mechanisms of deformation and failure in hcp Mg with a great accuracy [39,40]. Additionally, the NNp proposed in reference [40] detects the presence of the high pressure bcc-Mg phase with DFT accuracy [40]. The set of potentials studied in this work is listed in table 1. EAM potentials represent the most common model of atomic bonding in metals, in particular, with cubic symmetry. However, they do not work well for hcp metals and non-metallic systems, in which directional bonding is important, due to the approximation of the spherical shape of the electron cloud [52]. In the EAM-potential formalism, the total energy E EAM i of one atom i is given by the following expression [53]: where F is the embedding energy functional, ρ is the atomic electron density, φ is the pair potential interaction, and N is the number of neighbors inside a cutoff distance r c . Known the correct set of parameters for a given functional and density, EAM potentials have been proved to reproduce the structural properties of pure metals and metal alloys [53,54]. However, the functional form given by equation (1) can be modified to include additional contributions. After the development of the first Daw and Baskes EAM potentials described by equation (1) in the 80s, the TB-SMA models were the first modification of the EAM potentials applied to the study of pure Mg. The TB-SMA potentials have been proved to reproduce accurately elastic and transport properties in alloys [55]. In the present work, we analyze the performance of two TB-SMA models, which can be described by the Finnis-Sinclair formulation of the EAM potential [56], given by the following expression [50]: where A, p, r 0 , ξ and q are fitting parameters.
As we can see in figure 1, in the last 5 years, all efforts have been made in developing AD, MEAM and NN potentials. MEAM and AD potentials take into account the angular dependence of metallic bonds and are suitable for systems containing different types of bonding. Although these potentials are more computationally demanding than EAM potentials, they can improve the performance for alloys, especially with chemical species that adopt different lattice structures in their pure materials. MEAM potentials, which imply a slight modification in equation (1), are given by the following general expression [57]: In the EAM potential, the density ρ is given by a linear superposition of spherically averaged atomic electron densities, while in the MEAM potential the termρ includes angulardependent terms [57]. This modified potential was proved to work well with cubic materials and alloys [57], and constitutes a noticeable improvement for non-cubic materials. However, several modifications have been implemented in the formulation of Daw and Baskes of the EAM potentials since its origin in the 80s to capture more mechanisms at the atomic level and reproduce macroscopic properties with improved accuracy. AD potentials represent a generalization of EAM potentials by extending interactions to several additional coordination shells with angular-dependent terms. Their functional form is given by the following equation: [58]: where the dipole and quadruple distortions of the local atomic environment are introduced through the functions u and w, while s and t stand for the three Cartesian coordinates. While empirical potentials, such as EAM, MEAM and AD potentials, try to model interatomic interactions based on some physics behind, NN potentials are based on machine learning algorithms and their model parameters are chosen in such a way that they provide the atomic energy and the force on each atom as a function of descriptors of their atomic environments with near the same accuracy as that observed in DFT calculations. NN potentials have had a strong impact on the study of structures composed of a limited number of chemical species thanks to their small errors in the prediction of energies and forces with respect to DFT calculations, although at expenses of higher computational costs, with traditional NNps being up to 60 times slower than EAM potentials and 12 times slower than MEAM potentials [59]. Although these potentials are highly accurate, their fitting is still a challenging task, especially for multi-element compounds or alloys. In the NNp proposed by Stricker et al [39], which follows the Behler-Parrinello formulation [60], the local atomic energy E i of the atom i is described as: where j runs over all the neighbors atoms of the atom i within a cutoff sphere with radius r c , n is an integer, r e is the nearest neighbor distance, S i j is an angular screening term and α n is a model parameters. Thanks to the similarities with the MEAM formalism, it is possible to capture physical knowledge about metallic systems directly in the descriptor [40].
In this work, all the potentials listed above were tested by comparing their predictions for the cohesive curves, elastic constants, phonon frequencies and generalized stacking fault (GSF) energies of bcc Mg. The LAMMPS software package [61, 62] was used for energy calculations. The phonon dispersion curves were calculated with the Phonopy and PhonoLAMMPS codes [63].

Results and discussion
Unlike the other pure metals, Mg does not melt out from the bcc-phase at ambient pressure. At ambient pressure, the bcc structure is not stable in the pure material, regardless of the temperature. The hcp → bcc phase transition occurs near 52 GPa at 0 K, with decreasing transition pressure with increasing temperature [12]. Consequently, the interatomic potentials proposed for pure Mg need to capture additional atomistic mechanisms not present in other metals, and the parameterization of the interatomic potential, i.e. the choice of the appropriate model parameters in the selected mathematical function, becomes a more challenging task. The parameterization of an interatomic potential, as well as the choice of the functional form of interatomic potential, determines the prediction of the mechanical, thermodynamic, and transport properties of the final material. Usually, interatomic potentials are fitted to the bulk cohesive energy, equilibrium lattice parameter and elastic constants at 0 K [64]. However, the number of fitting parameters is limited by the functional form of the interatomic potential and often they cannot capture all atomistic mechanisms that give rise to material-specific properties. In general, obtaining accurate interatomic potentials that reproduce phase transitions is a challenging task, especially in those materials in which the equilibrium cohesive energies of the two phases are close by. Phase transitions depend on multiple factors, such as pressure, temperature or concentration, while interatomic potentials are often fitted to a reduced number of atomic configurations calculated with DFT at 0 K. Consequently, the capability of an interatomic potential to reproduce phase transitions in a wide range of concentrations, temperatures and pressures is considered as a benchmark when all these factors have not been taken into account in the fitting procedure.
At room temperature and ambient pressure, Mg presents the hcp structure and the hcp → bcc phase transition is observed at large pressure. Therefore, the bcc phase is not of great interest in the pure material and most interatomic potentials available in the literature were fitted to study the hcp structure. In this work, we compare the performance of fourteen interatomic potentials   [16]. available in the literature in the bcc structure. This comparison is focused on the study of cohesive curves and elastic constants, which characterize the mechanical properties of the material, and phonon dispersion curves, which provide information about the stability of the structure at finite temperature. Cohesive energy curves contain important information about energetics of atomic interactions in the crystalline structure as a function of the lattice volume, and describe how strong atoms are bound together in the lattice, also determining the mechanical strength of the structure. Thus, small cohesive energies correspond to weak atomic bonds, while large cohesive energies are observed in materials with strong bonds and good mechanical strength, such as insulators or semiconductors. In this work, these curves were extracted from energy calculations with LAMMPS of the perfect bcc-Mg lattice structure with varying lattice parameter. All potentials predict the equilibrium lattice parameter and cohesive energy of the pure material with small deviations. However, as we can see from figure 2, the cohesive energy obtained from each potential differs under lattice compression or expansion, and this results in different atomic behavior under external stresses. Additionally, the cohesive curves at small volumes differ too much from the ideal curve, which results in large errors in the transition pressures. The transition pressures were computed with the common-tangent method. In this method, the transition pressure is obtained from the slope of the common tangents of the cohesive curves of the bcc and hcp phases, as shown in figure 3. With the interatomic potentials tested in the present work, the transition pressure varies from 9 to 140 GPa, meaning large deviations from the expected value: 52 GPa. These discrepancies in the prediction of the cohesive energies at small volume and transition pressures may be caused by the small difference in the cohesive energies of the hcp and bcc phases. However, it is critically important to model the difference in these cohesive energies correctly for a reasonable description of the phase transition and the behaviour of the bcc phase and, since the bcc phase is not often considered in the fitting, the atomic mechanisms at small volume are not well captured by the existing interatomic potentials. The particular phase transition experimented by this material would require a careful fitting of the interatomic potential to include possible contributions from anharmonicity in the vibrational modes and deviations in the electronic structure.
Similarly, the elastic constants were calculated with simulation boxes with increasing volume to analyze their dependence on the lattice parameter and the behavior of the bcc phase under positive and negative stresses. The reference values obtained with DFT calculations were computed with the lattice parameter of the pure bcc phase, the equilibrium lattice parameter of bcc Mg in 5 nm/5 nm Mg/Nb multilayers and the lattice parameter of bcc Mg at 45 GPa [16]. While the errors in the elastic constants can be acceptable at large atomic volume, compared to the typical volumes studied in the hcp phase, the relative errors in these quantities increase with decreasing volume due to the lack of information about the atomic mechanisms occurring at these volumes captured by the interatomic potentials. Additionally, certain potentials with smaller relative error in C 11 , ADp 1 , TB-SMA 2 and the EAM potentials, do not satisfy the Born mechanical stability criterion [65], C 11 > C 12 . Born's criterion is found to provide a good prediction of loss of stability under tension. Therefore, given these instabilities, we can conclude that the bcc phase is not well reproduced and the study of these structures with these potentials would not be reliable. Figure 4 shows bulk and shear moduli versus lattice parameter with the equilibrium volume and under intermediate and large pressure, and compared to DFT calculations. This intermediate point corresponds to bcc structures that have been proved to be stable in nano-multilayers [16]. While the elastic properties converge with all potentials at negative pressure, they provide highly variable values at high pressure. Additionally, NNp 1 requires special attention due to its behavior under small volumes: decreasing bulk modulus and elastic constants with decreasing volume. This may be due to the different way how this type of potentials is fitted. While empirical potentials are fitted to reproduce macroscopic properties, NN potentials are fitted to reproduce the atomic energies and forces calculated from DFT calculations. It is likely that configurations with the bcc structure have not been used in the fitting procedure and therefore this wrong behaviour would be caused by an inaccurate extrapolation. This problem was solved in NNp 2 , which did include configurations with the bcc phase in the training set. According to table 2, the MEAM potentials are the potentials that best predict the mechanical properties of the bcc-Mg phase. This table shows the equilibrium lattice parameter, cohesive energy and elastic constants of bcc Mg as calculated by the different interatomic potentials analyzed in the present work, together with the hcp → bcc transition pressure. These Table 2. Equilibrium lattice parameter, a eq , in Å; cohesive energy, E c , in eV; elastic constants and hcp → bcc transition pressure, P hcp→bcc , in GPa, predicted by the interatomic potentials studied in the present work. The transition pressure was calculated with the common-tangent method from the cohesive curves. DFT values were taken from references [12,66]. The last column indicates the average relative error in the values of each potential with respect to DFT values. The potentials marked with an asterisk in the last column satisfy the Born condition, C 11 > C 12 , for the three lattice parameters.
a eq E c a = a eq a = 3.347 Å a = 3.043 Å P hcp→bcc η avg (%)   values are compared with DFT results obtained by Su et al [16]. The elastic constants were calculated for different lattice parameters. The better accuracy of the MEAM potentials may be caused by the limitations of the EAM potentials, such as bad performance in the presence of directional bondings and in phase transitions, which have been slightly corrected by the incorporation of modified versions that have been applied to hcp-Mg in the last five years. The main problem of the NN potentials is that they still do not incorporate enough information about the bcc phase. The mechanical properties of the material, characterized by the shear and bulk moduli, have a significant impact on dislocation behavior, interfacial behavior and dislocation-interface interactions [19]. We calculated the GSF energy for the (110) and (112) planes along different directions to study main dislocation slip systems in bcc Mg. The GSFE energies are calculated as a function of the unitary displacement vector and compared with DFT calculations [19]. For these calculations, once the two bcc crystals are created, one of the crystals is displaced a normalized distance u and the positions of the atoms are relaxed in the direction perpendicular to the interface. However, the relaxation of bcc structures is not stable with all potentials and the GSFE energies could only been calculated for a reduced set of potentials. Figure 5 shows the GSF energies calculated for the two different planes. The maximum stacking fault energies are obtained with the MEAM 3 potentials, while the minimum values are obtained with EAM 3 . This means that the development of deformation mechanisms is more favourable with the EAM 3 potential. However, none of the GSFE curves obtained with the different potentials matches the DFT values for all shear directions and planes and, in general, EAM potentials tend to underpredict surface energies [67], so their behaviour is not associated with the atomistic mechanisms occurring in bcc Mg atoms that they can capture, but to their functional form. Additionally, this behaviour reflects that these EAM potentials are not the best option for the study of surfaces or interfaces, such as in Mg/Nb multilayers. Figure 6 shows phonon dispersion curves along the symmetry directions for bcc Mg structure in fully relaxed (a) and under pressure (b) and (c) states, shown for different representative potentials. According to DFT calculations [18] shown in figure 6 as black squares, imaginary phonons (negative frequency) are found at the N point in fully relaxed bcc Mg. The presence of these imaginary frequencies indicates that this phase is not mechanically stable in its relaxed form and transforms to the stable structure: hcp. However, these curves become stable at reduced volume and the phonon frequencies increase with increasing pressure. This means that the bcc phase is mechanically stable at higher pressure and the stable phase at reduced volumes will be the phase with lower cohesive energy. The potentials studied in the present work do not reproduce the presence of imaginary phonons in the N -Γ direction; however, they are found in the H − P direction in the Brillouin zone with the EAM, ADP and NN potentials, and even with bcc Mg under high-pressure. These means that these constrained bcc phases cannot be reproduced with these potentials, as the hardening of the phonon branches is necessary to observe the stabilization of the bcc phase at high pressure or in regions close to Mg/Nb interfaces. On the other hand, the MEAM potential do not present imaginary phonons, which means that this structure is stable, although does not reproduce the atomic mechanisms present in Mg. However, empirical potentials are often known by their constraint on the prediction of both elastic and phonon properties with a great accuracy at the same time [67,68], and therefore the choice of the magnitudes used in the fitting procedure should be only the relevant ones for the phenomena of interest at expenses of a loss of accuracy on the other magnitudes.

Conclusions
The stabilization of metastable phases and understanding their properties at the atomic level is crucial to develop new materials with potential industrial applications. MD simulations are an ideal tool to study atomic-level mechanisms with a great accuracy. Unlike DFT calculations, in which the computer limitations do not allow for the study of large systems and long simulation times, MD simulations can be run with millions of atoms and time steps. However, these simulations require previous knowledge of interatomic potentials that reproduce the interatomic interactions. The development of these potentials is a challenging task, while also biased towards reproducing the properties of thermodynamically-stable phases.
Mg, as one of the most promising light-weight materials for industrial applications, adopts the hcp structure under normal conditions, resulting in low ductility. However, the stabilization of Mg in its bcc structure improves the mechanical properties of the hcp phase, but further studies on this phase are still necessary. The use of MD simulations would be a magnificent tool to analyze the stabilization of this phase, but most interatomic potentials available in the literature have been fitted to study the hcp phase and, as we showed in this paper, they do not reproduce the bcc phase as well as desired. According to table 2, the MEAM 3 and MEAM 4 potentials, which have a similar parameterization, are the potentials that best predict the mechanical properties of bcc Mg.
At ambient temperature and pressure, pure bcc Mg can exist as a metastable phase in sufficiently thin multilayer systems with a stiff metal, such as Nb. Mg/Nb multilayers have been proved to exhibit significantly higher yield stresses and strains to failure, compared to pure Mg [16]. However, the effect of the Mg/Nb interactions in the metastability of the bcc phase and the deformation mechanisms in the multilayers are still not well understood. Therefore, the development of new interatomic potentials for the Mg-Nb interaction are still necessary for the study of these multilayers in-depth at room temperature. Additionally, new advanced potentials, such as machine learning potentials, fitted with properties of the bcc structure are necessary to study bcc Mg at large pressure.