Atomistic Design of Favored Compositions for Synthesizing the Al-Ni-Y Metallic Glasses

For a ternary alloy system promising for obtaining the so-called bulk metallic glasses (BMGs), the first priority issue is to predict the favored compositions, which could then serve as guidance for the appropriate alloy design. Taking the Al-Ni-Y system as an example, here we show an atomistic approach, which is developed based on a recently constructed and proven realistic interatomic potential of the system. Applying the Al-Ni-Y potential, series simulations not only clarify the glass formation mechanism, but also predict in the composition triangle, a hexagonal region, in which a disordered state, i.e., the glassy phase, is favored energetically. The predicted region is defined as glass formation region (GFR) for the ternary alloy system. Moreover, the approach is able to calculate an amorphization driving force (ADF) for each possible glassy alloy located within the GFR. The calculations predict an optimized sub-region nearby a stoichiometry of Al80Ni5Y15, implying that the Al-Ni-Y metallic glasses designed in the sub-region could be the most stable. Interestingly, the atomistic predictions are supported by experimental results observed in the Al-Ni-Y system. In addition, structural origin underlying the stability of the Al-Ni-Y metallic glasses is also discussed in terms of a hybrid packing mode in the medium-range scale.

experimental results, however, their predictions have often shown some limitations. To the present authors' view, the limitations of the proposed criterion/rule are from their starting bases, e.g., the deep eutectic rule was based on the equilibrium phase diagram and the size difference rule was based on the atomic sizes of the constituent elements, and these starting bases have some restrictions in well reflecting the internal characteristics of the alloy system concerned 3 . The key is therefore to seek for a valid starting base and further develop an approach capable of clarifying the metallic glass formation in the specific alloy system. From a physical viewpoint, the interatomic potential of an alloy system is able to reasonably describe the major interactions involved in the system. Therefore, if a realistic interatomic potential is constructed and known, most of the physical and chemical properties of the system, including those related to the BMGs, can be deduced through relevant computations and simulations 3,12 . In the family of BMGs developed so far, the Al-based BMGs constitute a significant member, as they show unique properties in many aspects, such as high specific strength and even good ductility 13 . Among Al-based BMGs, Al-TM-RE (TM = Ni, Cu, Fe, Co, etc.; RE = Y, Ce, Gd, La, etc.) systems are found to be the most promising 14 . However, experimental observations of glass formation in the Al-TM-RE systems are often found to be in conflict with currently available rules or criteria and an in-depth understanding is demanded 13,15 . In the present study, the Al-Ni-Y ternary alloy system is selected as a representative of the Al-TM-RE systems for developing the atomistic approach. We propose to take a newly constructed Al-Ni-Y interatomic potential as the starting base together with a relevant simulation route to develop an atomistic approach which is capable of not only clarifying the formation mechanism, but also predicting a favored glass formation region (abbreviated as GFR hereafter) in its composition triangle and an amorphization driving force (abbreviated as ADF hereafter) for each possible glassy alloy located inside the predicted GFR. The predicted GFR indicates the energetically favored alloy compositions, which could serve as guidance for the composition design in synthesizing BMGs. The predicted ADF, related to the energy difference between the glassy phase and the crystalline solid solution counterpart, could give hint to the readiness of metallic glass formation for a specific glassy alloy located in the GFR, and may somehow be correlated with the technical defined GFA by either obtainable size or applied cooling speed 6 . It should be noted that both the predicted GFR and ADF are derived from the Al-Ni-Y potential of the system, reflecting mainly in an energetic or thermodynamic aspect. In practice, using different glass producing techniques, one could obtain different experimentally identified glass formation regions. The technical defined GFA could also be influenced by dynamic factors, such as the viscosity of the liquid melt, atomic diffusivity in the liquid, overheating during producing, etc. Consequently, the practically observed glass formation regions and the so-called GFA could exhibit some fluctuation from the predicted/revealed characteristics.
In addition, the properties of the condensed matters are believed to correlate with their atomic structures 16,17 . For the Al-TM-RE systems, it is intriguing that adding a few atomic percent of TM or RE as a third element could dramatically affect the glass formation and properties of the base binary alloy systems [13][14][15] . This sensitive alloying effect has been explained by some thermodynamic or kinetic arguments, including the equilibrium phase diagram, Gibbs free energy and fragility 18 . These arguments, however, fall within the macroscopic domain, and a microscopic picture at an atomic-scale is needed 19 . Consequently, resolving the atomic structure, monitoring the delicate modification of the characteristic short-and medium-range orders with TM or RE addition, could improve the understanding of such a microscopic picture and help elucidate its underlying structural origin. Although several atomic structural models have been proposed in the past decades [20][21][22][23] , owing to the complexity and diversity of the internal interactions in various alloy systems 3 , the atomic packing details in metallic glasses are still matters of intense debate. Note that the knowledge concerning the atomic structure of the metallic glasses is not only of scientific interest, but also, if well-clarified, would lead to an improved understanding and in turn help in controlling the properties of the metallic glasses 16,24 . Consequently, besides the discussion of the prediction of GFR and ADF from the Al-Ni-Y potential, we will also discuss how to apply the interatomic potential to characterize the atomic structure of the Al-Ni-Y metallic glasses via relevant computations and simulations 3   where E i is the total potential energy of atom i, φ and ψ are called here the pair term and density term, respectively. r m1 and r m2 are the knots, and r c1 and r c2 are the cutoff radii of the pair term and density term, respectively. n 1 and n 2 are indices that should not be less than 3 and 5, respectively, to avoid discontinuity of the high derivatives. A 1 , p 1 , A 1m , and p 1m and A 2 , p 2 , A 2m , and p 2m are another eight adjustable potential parameters. The potential parameters are determined by fitting to the physical properties such as cohesive energies, lattice constants, elastic constants and bulk moduli of the elements as well as of the stable or virtual intermetallic compounds in each of the binary systems. We now summarize the fitting results of the Al-Ni-Y potential. Table 1 displays the six sets of potential parameters for the Al-Ni-Y system. Table 2 gives the reproduced lattice constants, cohesion energies, elastic constants and bulk moduli for the stable structures of Al, Ni, and Y together with the corresponding experimental values [30][31][32] . It can be seen that the constructed potential can well reproduce the physical properties of the metals. To ensure that the potential could reflect the realistic atomic interactions involved in the systems, we have included a number of complicated compounds that are identified in the experiments (or found in the phase diagram) 33-39 , e.g., oP16-Al 3 Ni, hP5-Al 3 Ni 2 , cF24-Al 2 Y, tP20-Al 2 Y 3 , oP12-AlY 2 , cF24-Ni 2 Y, oP16-NiY 3 , and so on. Meanwhile, compounds with different compositions and crystallographic structures are also included, to reflect the atomic interactions in various chemical and structural environments. To aid the construction procedure, ab initio calculations are also conducted by the authors to derive the necessary physical properties of the compounds (detailed in Methods Section). virtual compounds in the Al-Ni, Al-Y 15 , and Ni-Y systems, together with those obtained by the ab initio calculations as well as some available experimental lattice constants. From Table 3, one can see that the physical properties reproduced by the interatomic potential match well with those derived by ab initio calculations or experiments, confirming that the constructed Al-Ni-Y interatomic potential can well describe the energetic and structural characteristics of the alloy phases in the system. Another evaluation test for the relevance of the potential is to determine whether the potential can describe the atomic interactions at non-equilibrium conditions. A common practice is to derive the equation of state (EOS) from the potential and then compare it with the frequently used EOS in this field, i.e., the Rose equation, which is considered to be universal for all categories of solids 40 . Figure 1 shows the EOSs derived from the constructed Al-Ni-Y potential and the corresponding Rose equations for fcc-Al, fcc-Ni, hcp-Y, cP4-AlNi 3 , cF24-Al 2 Y, and cP4-NiY 3 , respectively 15 . Note that the EOSs have not been used in the fitting process and can therefore be considered as an external measurement for the relevance. From Fig. 1, one can see that the EOSs derived from the constructed potential agree well with the corresponding Rose equations, suggesting that the potential is relevant in describing the atomic interactions, even if the system is far from the equilibrium state. Meanwhile, the total energies are smooth in the whole range, without any 'jumps' or discontinuities, thus avoiding the appearance of non-physical behaviors in the simulations.
Glass formation region of the Al-Ni-Y system. We now take the constructed Al-Ni-Y n-body potential as the starting base to develop an atomistic approach capable of clarifying the metallic glass formation.  To set up a relevant simulation route and model used in the simulations, here we first review the up-dated experimental results related to the metallic glass formation. Summarizing the experimental observations from various glass producing techniques, such as liquid melt quenching, ion beam mixing and solid-state amorphization, it is found that under these non-equilibrium glass producing techniques, the resultant alloy phase is either solid solution or glassy phase (i.e., in a disordered state), but not otherwise. This is because under these producing techniques, during which an effective cooling speed is estimated to be from 10 2 to 10 13 K/s, the available kinetic conditions are extremely limited and would certainly retard the ability of those complicated structured intermetallic compounds to nucleate and/or grow 3 . It follows that in the non-equilibrium producing process, the competing phase with the glassy phase is only the solid solution, which has one of the three simplest structures, i.e., fcc, hcp, or bcc, whereas the rather complicated intermetallic compounds, if exist, are excluded in this competition 3,41,42 . The above outcome from experimental observations suggests that predicting the GFR of a ternary alloy system can be converted into a scientific issue of splitting its composition triangle into two different types of regions, energetically favoring the glassy alloy and solid solution, respectively 3 . As shown above, a realistic Al-Ni-Y n-body potential is constructed and it governs the energetic states of all of the alloy phases, including the solid solution and glassy phase. It follows that the above scientific issue could be resolved by applying the constructed Al-Ni-Y potential to perform relevant atomistic simulations, in which solid solution models are used to compare the relative stability of the solid solution and its disordered counterpart over the entire composition triangle (detailed in Methods Section).
Accordingly, we have set up the Al x Ni y Y 1−x−y solid solution models over the entire composition triangle. The constructed models are then allowed to evolve under Monte Carlo (MC) simulations [43][44][45] within the isothermal-isobaric ensemble at 0 Pa and 300 K for 2 million steps, which is testified to be sufficient for the models to be fully relaxed (detailed in Methods Section). After sufficient simulation time, the initial solid solution models reach a relatively stable state, i.e., the drifts for all of the related dynamic variables are negligible. Inspecting the three-dimensional atomic configurations and pair correlation functions, it is revealed that when varying the compositions, the Al x Ni y Y 1−x−y solid solution models generally exhibit two different states, either preserving the initial crystalline state, or collapsing into a disordered state, corresponding to the formation of a glassy alloy. Simulations over the entire composition triangle allow locating the GFR of the Al-Ni-Y system, shown in Fig. 2. One can see that the composition triangle is split into two types of regions by three critical solubility lines AB, CD, and EF. When an alloy composition lies beyond AB, CD, and EF and moves towards one of three corners, the crystalline solid solution could remain stable. These regions are thus classified as the crystalline regions. When the composition falls in the central hexagonal region enclosed by ABCDEF, the solid solution becomes unstable and would spontaneously collapse into the disordered state. This lattice collapse, or say the crystalline-to-amorphous transition, is a result of the relaxation of the atomic-level stress when an adequate amount of solute atoms dissolve into the solvent lattice. Once the solute concentration exceeds a critical value, the severe lattice distortion could trigger a collective collapse of the crystalline lattice, turning into a more topologically stabilized disordered state than the crystalline state. This hexagonal region enclosed by ABCDEF is defined as the GFR of the Al-Ni-Y system. Within the GFR, the formation of the Al-Ni-Y metallic glasses is energetically favored. In addition, in the vicinity of the boundary of GFR, there are several models that reside in an ordered-disordered coexisting state. The formation of such coexisting state could be attributed to the stability of the solid solution and glassy phase being relatively close, leading the ordered and disordered states to coexist. Figure 2 also indicates that for the Al-Ni and Al-Y binary sub-systems, the addition of a third element of Y or Ni helps the metallic glass formation. In the Al-Ni sub-system, the glassy phase could be obtained within a composition range of 10-60 at% Ni, and a few percent addition of Y would extend the GFR from one end to another in the composition triangle, suggesting that the ternary Al-Ni-Y metallic glasses could be formed at any combinations of Al and Ni. A similar case is also observed in the Al-Y sub-system, indicating that adding a third element of Ni also helps the metallic glass formation. Apparently, such a sensitive alloying effect obtained by adding a minor third element is of practical importance in synthesizing metallic glasses.
Composition optimization for glass formation. We now further discuss the composition optimization for glass formation in the Al-Ni-Y system by calculating an amorphization driving force (ADF) for each glassy alloy inside the predicted GFR.
In practice, it has been observed that the Al-enriched metallic glasses with Al content over 60 at% feature good ductility 13,14 ; we therefore focus on the Al-enriched corner, i.e., Al x Ni y Y 1−x−y (x = 60-100 at%, y = 0-40 at%) alloys for composition optimization. From an energetic point of view, the formation enthalpy difference between the glassy and solid solution phases (Δ H glassy − Δ H s.s. ) could serve as an amorphization driving force (ADF) for the specific alloy, whereas the formation enthalpy of the solid solution phase (Δ H s.s. ) could act as a resistance against amorphization. We propose to define a parameter, γ, namely the normalized ADF for a specific glassy alloy inside the predicted GFR, and it can be expressed by  For Δ H s.s. , all of the solid solution models need to be relaxed to reach their respective minimum energies at various compositions. In relaxation, the MC box is set to be fixed, i.e., retaining the initial crystalline structures, whereas only atom displacement takes place in the process. The models are then relaxed for an adequate time period. Assuming that E min is the minimum energy per atom of the Al x Ni y Y 1−x−y solid solution, the formation enthalpy Δ H s.s. of the solid solution could then be expressed by The contour map of γ for the Al x Ni y Y 100−x−y (x = 60-100 at%, y = 0-40 at %) alloys is plotted in Fig. 3. Inspecting Fig. 3, it can be seen that when the compositions are located inside the area marked with red dots, their corresponding γ parameters are larger than those residing outside. From further calculations, a maximum value of γ can be deduced at a stoichiometry of Al 80 Ni 5 Y 15 , and the surrounding sub-region marked with deep red dots could be considered as the optimized compositions for the Al-Ni-Y metallic glass formation. This means that from an energetic point of view, if an Al-Ni-Y alloy is designed with its composition located in this small sub-region, the obtained metallic glass would probably be prominently stable, and that from a kinetic point of view, this glassy alloy might be more easily obtained in practice than those located outside the sub-region.

Discussion
Comparison with the experimental observations. We now discuss the comparison of the predictions for the GFR and ADF by the atomistic approach with the experimental observations reported so far in the literature.
First, to validate the predicted GFR for the Al-Ni-Y alloy system, we have collected the available data 13,14,46 and marked them in Fig. 2 46 . One can see from Fig. 2 that the compositions of these experimentally obtained metallic glasses mostly fall within the predicted GFR, indicating that the proposed atomistic approach is quite reasonable in predicting/ determining the GFR for the Al-Ni-Y system, allowing one to conveniently design appropriate compositions for synthesizing the Al-Ni-Y metallic glasses. Meanwhile, as illustrated by the yellow shaded area in Fig. 2, the compositions of the experimentally obtained metallic glasses mostly fall within the Al-enriched corner, suggesting that this area is indeed promising for searching for the optimized compositions for forming metallic glasses.
Second, it is of interest to determine if the atomistic prediction of the ADF for individual Al-Ni-Y alloys could somehow be correlated with the technically defined GFA, which is frequently cited in the literature. Inspecting Figs 2 and 3, one can see that experimentally measured compositions of the obtained Al-enriched metallic glasses are densely distributed in the red area in Fig. 3, within which the normalized ADF (γ) parameters are calculated to be larger than those compositions sitting outside. Furthermore, it is found that the Al contents of the best glass formers determined in the experiments are typically around 85-86 at% [46][47][48][49] . Interestingly, these experimental compositions mostly fall within the deep red sub-region (nearby the stoichiometry of Al 80 Ni 5 Y 15 ), which is deduced to include the optimized compositions for the Al-Ni-Y metallic glass formation. It is noted that the pinpointed optimized compositions are defined by a small sub-region, but not a single number. This is because the above calculations/simulations are based on the constructed Al-Ni-Y n-body potential, in which a total of 6 × 15 parameters are fitted with the data from experiments or ab initio calculations, and the precision of the calculation/simulation results should be discussed with consideration of the inevitable errors involved. It is commonly accepted that in the atomistic simulations, the error is approximately 3-5%. A simple calculation then identifies a small composition sub-region around the maximum value of Al 80 Ni 5 Y 15 , and the compositions inside this small sub-region could be considered as the optimized compositions for the Al-Ni-Y metallic glass formation. From the above analyses, it can be seen that the experimental results are in support, in a qualitative or at least semi-quantitative manner, of the ADF predicted from the present atomistic approach for the Al-Ni-Y system. Moreover, as mentioned above, the predicted ADF represents the important energetic factor governing the metallic glass formation, yet it is not an exclusive one. For the technical defined GFA, it would be affected by a variety of factors in practice, and in addition to the energetic aspect, the liquid structure, viscosity, atomic diffusivity, as well as possible defects, impurities and microheterogeneities, etc., should all be considered. Therefore, the agreement in the present study between the atomistic prediction and technically measured parameters is quite reasonable. In other words, the predicted ADF provides a comparative measure, from an energetic point of view, for the GFA defined by some technical parameters, e.g., critical size or cooling speed.
In addition, it has been reported that the thermodynamic calculation method has also been employed to predict the glass formation in the Al-Ni-Y system 50 . However, it was found that the predictions exhibit considerable discrepancies with the experimental observations, especially in the Al-enriched corner. This might be attributed to the semi-empirical and phenomenological nature of the thermodynamic calculation method. The deviation between the thermodynamic predictions and experimental observations in the Al-enriched corner has also been revealed for other Al-TM-RE systems, such as the Al-Ni-(La,Ce) systems 50 . Moreover, a simple topologically-based prediction scheme has also been proposed to pinpoint the composition for Al-Ni-Y system and the predictions are quite acceptable 47 . However, this scheme is still based on a not so explicit starting base, i.e., cluster stability, and the predicted (local) optimal compositions are insensitive to the types of constituent atoms and similar for different alloy systems. Consequently, we have shown that the present atomistic approach, based on the interatomic potential, could be an improved approach in dealing with the issues of metallic glass formation in the ternary alloy systems.
Summarizing the above analyses, it is revealed that the atomistic prediction of the GFR and ADF in the Al-Ni-Y ternary alloy system is relevant and supported by the experimental observations. As the atomistic prediction is directly derived from the constructed Al-Ni-Y n-body potential, the relevant prediction could, in turn, provide additional evidence to the relevance of the constructed potential. Atomic Structure and stability of the metallic glasses. We then scrutinize the structural characteristics of the Al 80 Ni 5 Y 15 metallic glass, which features the largest γ, implying that it may have the largest phase stability as well.
The short-range order (SRO) is characterized first. The Voronoi tessellation method is employed, and cell faces with smaller than 5% of the average face area are excluded to minimize the degeneracy problem and thermal vibration effects 45 . As presented in Fig. 4a, it can be seen that the CN spectra in Al 80 Ni 5 Y 15 is distributed over a wide range, with three types of atoms covering different scopes of whole landscape, i.e., Ni dominates in CNs of 9-11 (average 9.8), Al dominates in 12-15 (average 13.3), and Y dominates in 16-19 (average 17.2). As the Al-Ni-Y system is a typical system comprised of small-, medium-and large-sized atoms, i.e., the Goldschmidt radii of Ni, Al and Y are 0.124 nm, 0.143 nm and 0.180 nm, respectively, thus, naturally, the larger radii of Y permit accommodate more atoms in the neighboring shells and lead to larger CNs, followed by Al, and then Ni. When forming metallic glasses, clusters of different sizes can better coordinate in space and efficiently fill the sites in the disordered structure, leading to the increase in packing density as well as the enhancement of phase stability. Furthermore, inspecting the distribution of Voronoi clusters, it is found that Ni and Y, as solutes, are mainly, but not strictly, surrounded by solvent Al as nearest-neighbors and mostly form Kasper-type clusters, such as Ni-centered < 0, 2, 8, 0> and Y-centered < 0, 1, 10, 6> . These solute-centered clusters can be considered as the building blocks in Al 80 Ni 5 Y 15 . However, this is comparable but slightly different from the scheme of solute-centered quasi-equivalent packing 21 , as no particular "solute-solute avoidance" is observed in Al 80 Ni 5 Y 15 . In previous studies, for the Al-Ni-La metallic glasses, such as Al 89 Ni 5 La 6 with a larger Al content of 89 at%, the solute Ni and La are almost totally avoided, exhibiting a large extent of chemical SRO 51 . However, in the present study, Ni and Y atoms are more close to random packing, with a Warren-Cowley parameter ~0. The packing of Ni and Y is revealed in Fig. 4b, where the neighboring linkages are displayed. This behavior can be attributed to the fact that the mixing enthalpy of Ni-Y is also largely negative, i.e., − 31 kJ/mol 11 , which is comparable to those of Al-Ni and Al-Y, i.e., − 22 and − 38 kJ/mol, respectively 11 . Therefore, Ni and Y are not necessarily neighbored exclusively by Al and Ni-Y neighboring is also stabilized energetically. Meanwhile, as the total content of Ni and Y are up to 20 at% in Al 80 Ni 5 Y 15 , which is larger than typical solute-lean metallic glasses 51 , a certain extent of neighboring between Ni-Y is also expected from the nominal composition.
Characterizing SRO alone, neglecting the heterogeneous correlations among various types of clusters that extend to the next level of hierarchy, i.e., medium-range and beyond, is insufficient to explain the structure-property relationship of metallic glasses 52,53 . A microscopic picture of medium-range order (MRO) can then be resolved by the percolated network formed by the solute-centered clusters. In MRO, the clusters adopt dense arrangements, by vertex-(VS), edge-(ES), face-(FS) and tetrahedra-sharing (TS) linkages 54 . In Fig. 5a,b, typical packing of neighboring clusters around Ni and Y clusters is illustrated, indicating that the VS, ES, FS and TS linkages are collectively hybridized in achieving the efficient packing in the metallic glass. In Al 80 Ni 5 Y 15 , the Ni-centered clusters typically have 10-13 solute-centered clusters as neighbors, whereas the Y-centered clusters often have 15-18 clusters around, suggesting more compact packing around Y clusters. Scrutinizing the packing of solute-centered clusters, a hybridized packing mode is observed, of which the icosahedral-like (five-fold) and fcc-like (six-fold) arrangements are most prevalent, as highlighted by the black dotted lines in Fig. 5a,b, and the higher-coordinated fcc-like arrangements are more populated around Y clusters. This observation agrees with our previous findings for the other Al-TM-RE systems 15 . The hybridized packing mode cannot be fully covered by previously proposed structural models, which often suggest a single packing mode for simplicity 16,20,21 . Our findings suggest that despite these structural models capturing the efficient packing nature of MRO, they are still providing an idealized, or over-simplified, picture. For real-world metallic glasses, due to the distinct atomic sizes and chemical interactions involved, MRO has its intrinsic complexity, just as revealed by the hybridized packing mode for the Al-Ni-Y metallic glasses. A close-up view of the highlighted icosahedral-like and fcc-like packing is further presented in Fig. 5c,d, respectively. This complexity of MRO creates more opportunities for space tiling, facilitates constituting a well-percolated network that can serve as the reinforced 'backbone' , and eventually leads to the improved stability of the metallic glasses 49 .
To summarize, we have shown that by taking a proven realistic Al-Ni-Y n-body potential as the starting base to conduct series simulation under a relevant route over the entire composition triangle, an atomistic approach is developed. The atomistic approach clarifies that the physical origin of the metallic glass formation is the collapsing of solid solution when the solute concentration exceeds the critical value and predicts an energetically favored glass formation region (abbreviated as GFR) of the ternary alloy system. Moreover, the atomistic approach further predicts an amorphization driving force (abbreviated as ADF) for each alloy located within the GFR. The predicted ADF could be a comparative measure, from an energetic perspective, of the so-called glass-forming ability (known as GFA), which has long been used in the field, yet was defined by some technical parameters, e.g., critical size of the obtainable glass. In addition, structural analysis reveals that a hybridized packing mode exists in the Al-Ni-Y metallic glasses, manifesting the intrinsic complexity of MRO and relevantly interpreting the phase stability of the Al-Ni-Y metallic glasses. Obviously, the approach presented here could serve as an important guide in designing appropriate alloy compositions for synthesizing ternary BMGs and in finding new possible BMG formers as well.

Methods
Construction procedure of the Al-Ni-Y potential. In the present study, a set of the Al-Ni-Y n-body potentials is constructed under the smoothed and long-range TB-SMA formulism. Concerning the atomic interactions in the Al-Ni-Y system, there should be six sets of potentials, i.e., three sets for Al-Al, Ni-Ni and Y-Y, and three sets for Al-Ni, Al-Y and Ni-Y. The latter three sets are referred to as cross potentials, which describe interactions between dissimilar atoms. For each set of potential, there are 15 potential parameters to be fitted, i.e., < p 1 , A 1 , r m1 , n 1 , p 1m , A 1m , r c1 , p 2 , A 2 , r m2 , n 2 , p 2m , A 2m , r c2 , r 0 > . Specifically, the potential parameters of Al-Al, Ni-Ni and Y-Y are determined by fitting them to the experimental properties of the Al, Ni and Y metals at 0 K, such as the cohesive energy, lattice constants, elastic constants and bulk moduli [30][31][32] . The parameters of Al-Ni, Al-Y and Ni-Y cross potentials are determined by fitting them to the physical properties of several stable or virtual intermetallic compounds in each of the binary systems. However, there are always few available property data of the related compounds that could be used in the fitting procedure. To solve this problem, ab initio calculations are then conducted by the authors to derive the necessary physical properties of the compounds and assist the potential construction.
The ab initio calculations are carried out using the Cambridge serial total energy package (CASTEP) 55,56 based on density functional theory (DFT). During the ab initio calculations, the exchange-correlation energy functional is described by the established Perdew-Wang (PW91) version of the generalized gradient approximation (GGA) 57,58 , and the ion-electron interactions are treated by the ultrasoft Vanderbilt pseudopotential scheme (US-PP) 59 . The cutoff energy is chosen to be 500.0eV, and the Brillouin zone is Scientific RepoRts | 5:16218 | DOI: 10.1038/srep16218 sampled using the Monkhorst-Pack method 60 with nearly constant k-point densities for each calculation, roughly equivalent to a 12 × 12 × 12 mesh for a conventional fcc unit cell. These parameters are shown to be sufficient for convergence of the calculations for the structures in this work. Each structure is firstly optimized with respect to the external degree(s) of freedom as well as the internal degree(s) of freedom (if any) of the unit cell as permitted by the space-group symmetry of the crystalline structure. The lattice constants and total energies of the optimized structures are then obtained, and the elastic constants and bulk moduli can be derived by fitting the stress-strain relationship.
Simulation model. To compare the relative stability of the solid solution and its amorphous counterpart, the solid solution models are commonly employed as the simulation models 3 . Two types of solid solution models, i.e., fcc and hcp, are established for the Al-Ni-Y system. For the fcc models, the [100], [010], and [001] crystalline directions are parallel to the x, y, and z axes, respectively. For the hcp models, the [100], [120], and [001] crystalline directions are parallel to the x, y, and z axes, respectively, and the hcp models can be considered to be built of equivalent orthorhombic cells. Each solid solution model consists of 2916 atoms, with periodic boundary conditions adopted in three Cartesian directions. To conduct a thorough investigation on the entire compositional phase-space, we have established Al x Ni y Y 1−x−y models over the entire composition triangle. The lattice parameters of the Al x Ni y Y 1−x−y solid solution models start out following the Vegard's law. In setting up the solid solution models, solute atoms are added by random substitution of a certain number of solvent atoms to reach a desired concentration.
Monte Carlo simulation. The established solid solution models are then allowed to evolve under Monte Carlo (MC) simulations within the isothermal-isobaric ensemble at 0 Pa and 300 K for 2 million steps, which is shown to be sufficient for the models to be fully relaxed. In order to minimize the possible random error, the MC simulations have been independently performed for five times and the calculation results are averaged. During the simulations, structural changes occurring in the models are monitored by the three-dimensional (3-D) atomic configuration, which can visually reflect the state of the system, and the pair correlation function, which is commonly recognized as firm evidence to identify the glassy phase.
The details of the present MC simulations are briefly described as follows. In MC simulations, there are two types of moves: atom displacement and box deformation. For atom displacement, the simulation system can be treated as a canonical ensemble at constant NVT, i.e., the shape and volume of the box are fixed. An atom i is chosen at random and given a uniform random displacement, where Δ E = E trail − E is the energy change resulting from the trial move of atom i. k B and T are the Boltzmann constant and the temperature of the simulation system, respectively. For box deformation, the fractional coordinates of atoms in the box are fixed. The state of the simulation box is defined by the matrix h, from which the trial state h trail is generated by the following relations, where Δ E = E trail − E is the energy change that results from the trial move of the box. σ is the applied stress tensor. Ω trail is the volume of the box in the trial state. When only hydrostatic pressure is applied, P move One can see from Eq. (9) that there are 6 independent elements in the strain tensor. To monitor the internal stress or pressure of the system, the 6 strains are separately applied to the box in a MC step and therefore the internal stress can be synchronously computed by For a system consisting of N atoms, there are a total of N + 6 trial moves; N moves for atoms and 6 moves for the box. The N + 6 trial moves are randomly chosen to be carried out in an MC step. In addition, if a move is rejected, the old state is recounted. The amplitudes of trial moves δ A and δ B are also adjusted so that the acceptance ratios of the trial moves are kept around 50%.