The Purpose of the Derivation Below Is to Show How the Resolved Shear Stress (And Yield

A. Development of the Cu-Zr-Al ternary EAM potentials

As stated in the manuscript, the interatomic potentials employing the EAM formalism [1] were obtained by fitting the potential energy surfaces (PES) [2] of the Cu-Zr-Al system. According to Ref. [1], the potential energy of an embedded atom i can be expressed as: , where α and β are the element types of the center (embedded) atom i and the neighboring atom j. ρβ is the charge density function of β; Fα is the embedding function ofα; Φαβ is the pair potential function between α and β. Note that the charge density here is the “effective” density for EAM, not the true charge density. Total potential energy of the system is the summation of the potential energy of each atom. The charge density function and embedding function of each element, as well as the pair potential function for each element pair, are thus the target functions of potential fitting. The fitted functions of our EAM potentials are given in a tabulated/numerical form, withthe profile of the functions shown in Fig. S1.

Over 600 configurations (typically 96~128 atoms obtained via ab initio MD [3] at different temperatures) of Cu, Zr, Al, Cu-Zr, Cu-Al, Zr-Al, Cu-Zr-Al are employed in the potentialfitting, including crystal structures and intermetallics, crystals with interfacial and point defects, solid solutions, liquid phases, MGs quenched at different cooling rates, and liquid inherent structures, etc. Experimental elastic constants and phonon frequencies of selected crystals are also included in fitting. For each selected atomic configuration, high-precision ab initio calculations were performed to obtain the PES, where the Brillouin zone was sampled with 2×2×2 Monkhorst-Pack k-point grids. The modifiedpotfit program [4] was used for the EAM parameter fitting. The generated EAM parameters were further refined using a recursive method: (1) Classical MD using the existing EAM parameters were used to generate a set of atomic configurations. (2) Ab initio calculations were performed on these configurations, and the data were incorporated into the datasets describing the PES (energies, forces and stress tensors). (3) The updated datasets (energy landscapes) were refitted to obtain new, improved EAM parameters. The process was repeated until self-consistent results were reached. The EAM parameters were validated against a large set of experimental and ab initio data, such as lattice parameters, cohesive energies, enthalpies of mixing, elastic constants, etc.In the following section we will provide detailed information about the validation and performance of the as-developed EAM potentials.

Figure S1. As-developed EAM functions. (a) Charge density functions. (b) Embedding functions. (c) Pair potential functions.

B. Validation and performance of the EAM potentials

B1. Valicaiton of crystal structures: EAM, ab initio and experiments

First of all, we consider crystalline phases consisting of one or more of the three elements (Cu, Zr, Al). In Fig. S2, various intermetallic compunds in the Cu-Zr-Al system are labeled in the ternary phase diagram. These compounds, as potential competing crystalline phases during glass formation, are of special importance.Therefore, accurate description of these compounds is the first requirement for the EAM potentials.To achieve this goal, equations of state for the pure elements (Zr, Al, Cu) and intermetallic compounds of Cu-Zr, Cu-Al, Zr-Al, and Cu-Zr-Al were first obtained using ab initiomethods. The cohesive energies of the corresponding atomic configurations were then calculated and fitted using the EAM formalism.The fitting results are shown in Fig. S3 to Fig. S7.The comparison shows that for the sampled configurations along the equations of state, the developed EAM potentials can reproduce the cohesive energies quite well. The average difference of the cohesive energybetween EAM and ab inito calculations is less than 0.015 eV/atom.

Figure S2. A schematic phase diagram of Cu-Zr-Al ternary system, with known intermetallic compounds labeled byopen circles.

Figure S3. Validation of EAM potentials for pure elements, including fcc Cu, fcc Al, and three Zr crystalline phases. Lines represent calculations with our EAM potentials, and symbols are obtained from ab initio calculationsusingVASP (the same for Figs. S4 to S7).

Figure S4. Validation of EAM potentials for various Cu-Zr compounds.

Figure S5. Validation of EAM potentials for various Cu-Al compounds.

Figure S6. Validation of EAM potentials for various Zr-Al compounds.

Figure S7. Validation of EAM potentials for two Zr-Cu-Al compounds.

To further validate the as-developed EAM potentials, we proceed to solve crystal structures employing the energy minimization techniquewith the EAM potentials, The EAM-derived structures were then validated against either ab initio calculations or experimental data available in literature. Various properties including lattice parameters, cohesive energies (after minimization), elastic constants, etc. have been calculated and compared. It is worth noting that most experimental values are measured at room temperature, while the simulation results are for 0 K. This temperature effect may contribute to the differnces in the listed properties, but the effect should be very limited for the properties we compared below. For example, due to thermal expansion, the lattice constants at 0 K and 300 K are estimated to differ by only ~0.5%. The detailed comparisons are compiled in the following tables.

TableS1. Comparison of ground state properties of elemental Zr, Cu, and Al.

a B. Olinger and J.C. Jamieson, High Temp. High Pressures 5, 123 (1973)

b J.R. Morris, Y.Y. Ye, K.M. Ho, et al., Phil. Mag. A 72, 751 (1995)

c A. Serra and D.J. Bacon, Phil. Mag. A 73, 333 (1996)

d Present work, ab initio calculation.

e C. Kittel, Introduction to Solid State Physics (whiley-Interscinece, New York, 1986)

f Metal Reference Book, 5th ed., edited by C.J. Smith (Butter-worth, Longdon, 1976)

g G. Simons and H. Wang, Single Crystal Elastic Constants and Calculated Aggregated Properties (MIT press, Cambridge, MA, 1977)

h R.W. Siegel, J. Nucl. Mater, 69&70, 117 (1978)

i H.E. Schaefer, R. Gugelmeiser, M. Schmolz, and A. Seeger, Mater. Sci. Forum 15-18, 111 (1987)

TableS2. Comparison of cohesive energies for other Zr phases.

TableS3. Comparison of cohesive energies for other Cu phases.

TableS4. Comparison of cohesive energies for other Al phases.

Tables S1 to S4 demonstrate that our EAM potentials can describe the properties of the three individual elements quite well, including the ground state properties, as well as the order of the various crystalline phases in terms of cohesive energy. The cohesive energy of the ω-Zr appears slightly lower than the equilibrium α-Zr (hcp), but similar results using ab initio calculations have been reported in several papers (see Ref. [5] and the references therein), and it has been suggested that ω-Zr could be the ground state of Zr at zero temperature and pressure, while α-Zr is the ground state at room temperature and ambient pressure in air. This is consistent with our simulation results. In addition to elemental properties, the properties of the intermetallic compounds predicted by EAM are summarized in Table S5.

TableS5. Comparison of compound properties.

* Lattice parameters of intermetallic compounds are found from Inorganic Crystal Structure Database (ICSD, e.g.,
** Note that thesephases are hypothetical and may not exist in experiments (phase diagrams).These hypotehtical structures are constructed mainly for potential fitting purposes. Lattice parameters of these structures are obtained through ab inito calculations.

It is seen that the match between EAM,ab initiocalculations, and experimentsin Table S5 is satisfactory. Moreover, the consistency of cohesive energies (in Tables S1 to S5) for both the elements and the compounds indicates a good match of formation enthalpy. In fact, the formation enthalpies can also be compared with experimental values in literature [6]. For example, using simulation (EAM or VASP), the formation enthalpy of CuZr2 (I4/mmm) is determined to be -14.2 kJ/mol, consistent with experimental measurements: -11 kJ/mol in Ref. [7] and -17 kJ/mol in Ref. [8]. For ZrAl (Cmcm), the calculated value is -43.4 kJ/mol, while the experimental value is -45 kJ/mol [9]. For CuAl (I2/m), the comparison is: -23.2 kJ/mol (EAM), -22.2 kJ/mol (VASP), and -20 kJ/mol (experiment, [10]). These examples demonstrate that not only the elemental cohesive energies, but also the chemical affinity between the constituent elements can be faithfully represented by the EAM potentials.

It would be interesting to compare the full icosahedral clusters in our simulated Cu46Zr54 MG with the local atomic packing in CuZr2 crystals, as they are similar in composition.Two artificial structuresof CuZr2 (with prototypes of NiTi2 and Al2Cu, respectively) are also considered [11], in addition to the stable CuZr2 with prototype MoSi2(see Table S5). The two artificial structures (including their internal coordinates and lattice parameters) of CuZr2 are obtained by direct energy minimization of the prototype configurations, and the lattice constants compare well with those reported in Ref. [11].As discussed in the paper, the full icosahedral cluster in Cu46Zr54 MG has a complete five-fold environment, but the shape is not a perfect icosahedron. Moreover, the chemical make-up and arrangement of the shell atoms also vary from site to site (Fig. 2c in the paper). The average Cu-Cu bond length is d~2.54Å, while the average Cu-Zr bond length is d~2.85 Å (Fig. 1 in the paper). Inthe CuZr2(MoSi2), in comparison,the center Cu is surrounded by 8 Zr atoms at d=2.85 Å, 4 Cu atoms further away at d=3.22 Å (as opposed to the nearest-neighbor Cu-Cu pair with d=2.54 Å in the MG),and another 2 Zr atoms at d=3.88 Å. More importantly, the bond symmetry is predominantly four-fold, not five-fold. In the CuZr2 (Al2Cu), Cu atoms are surrounded by 2 Cu and 8 Zr atoms, withCN=10 and the Voronoi index of <0,2,8,0>. Such a local environment is apparently different from the full icosahedra in MGs. The CuZr2(NiTi2) is interesting, becausesimilar to the full icosahedra in the MGs, the Cuin this crystal is also surrounded by 12 atoms, forming icosahedral local packing with complete five-fold environment. However, a close inspection reveals that the local atomic packing of this structure is still different from the full icosahedral packing inMGs in several aspects. First, the Cu-Cu bond length in CuZr2 (NiTi2) is d=3.1 Å, much larger than the Cu-Cu bond length in MGs,i.e. Cu and Cu are not in direct “contact”in CuZr2 (NiTi2). Second,in CuZr2 (NiTi2), the composition in the shell is always 3 Cu and 9 Zr atoms; while in the MG, 4 Cu 8Zr and 5 Cu 7 Zr are more favorable (see Fig. 2c in the paper).Third,in CuZr2 (NiTi2), the 3 Cu atoms in the shell and the center Cu atom always form a tetrahedron; while in the MG, Cu atoms in the shell can have various distributions.Fourth,in CuZr2 (NiTi2), Zr atoms are also surrounded by 12 atoms with icosahedral local packing;while in MG, the average CN around Zr is larger than 13.On top of all the above, one should note a general difference between the MG and crystals: the arrangement of shell atoms is ordered and identical in crystals, but disordered and diverse in MGs.

Aside from geometric considerations, one can also comparethe cohesive energy and stability of these competing compounds. A direct observation from Table S5 is that the three crystalline phases of CuZr2 have correct order and similar energy differences in both ab initio and EAM calculations. The MoSi2 structure has the lowest energy and is the stable phase of CuZr2 in experiments. This is consistent with the ab initio calculations in Ref. [11]. Although the artificially constructed CuZr2(Al2Cu and NiTi2) can survive direct energy minimization in simulations, they are not stable upon heating. For example, by increasing temperature to 500 K (below Tg), the CuZr2 (NiTi2) deconstructs rapidly in simulation, i.e.it fails torecover to the initial NiTi2structureupon subsequent energy minimization.

Based on the above discussion, we conclude that our simulated structure is truly amorphous;the full icosahedral cluster (with its fraction increasing with decreasing cooling rate) is neither a crystal nor a sub-critical nucleus of possible competing crystals.The potentials did not unphysically favor this local packing.

In addition to the low-energy crystalline states compared and studied above, many other configurations are also included in potential fitting, such ascrystals with interfacial and point defects, solid solutions, liquid phases, MGs (including Cu-Zr binary and Cu-Zr-Al ternary) quenched at different cooling rates, and liquid inherent structures, etc.This is important because fitting involving highly activated configurations (e.g. liquids and inherent structures) makes the EAM potentials more robust and transferable, especially for studying metastable systems such as MGs and supercooled liquids. In othe words, the potentials should be able to describe not only the PES around the crystalline phases, but also the regions around the activated states, with energies and forces fall within a reasonable error bar. The fitting results of all these(~600) configurations are summarized in Fig. S8 below.

Figure S8. Validation of the EAM potentials: comparison of pressure and energy with those calculated usingab initio simulations.

Fig. S8 demonstrates that the cohesive energies calculated by EAM have a high accuracy, with an average error of ~0.015 eV/atom compared with ab initio calculations. Note that the configurations with relatively large error (e.g. the largest 0.05 eV/atom) are mostly highly activated states such as liquids, inherent structures, or structures under high pressures. As we have shown before in Figs. S3 to S7 and Tables S1 to S5, the average error is much smaller for crystalline phases (typically0.015 eV/atom).The satisfactory match of cohesive energies for various liquid states also indicates that the heat of mixing (of high temperature liquids) compares well between the EAM and theab initio results. Here we note that the predicted values of heat of mixing using Miedema’s model (Cu-Zr: -23 kJ/mol, Zr-Al: -44 kJ/mol, Cu-Al: -1 kJ/mol [6]) have been widely used to evaluate the chemical affinity in Cu-Zr-Al alloy. However, the model prediction has limitations, especially for Cu-Al pair, where the electronic interaction could be unusual, leading to bond shortening (see Supplemental Material II). Therefore, in certain cases, the chemical affinity of Cu-Al could be underestimated by the model, thus accurate evaluation using ab initio calculation (or experiment [10]) is necessary.

The pressure error (ΔP) indicatesthe difference between the calculated pressuresofa configuration using EAM and ab initiosimulations respectively.Since the typical bulk modulus for this system is on the order of 100 GPa, the pressure error(ΔP) of 0.35 GPa can be roughly converted to a density error(Δρ) of 1%.

B2. Validation of liquid/glass structures

B2.1 Comparison of pair distribution functions: EAM vs. ab initio

We also compared the structures of MGs (as well as their high-temperature liquds) simulated using EAM MD andab inito MD. The partial pair distribution functions (PDFs) are shown for Cu46Zr54in Fig. S9, and for Cu46Zr47Al7in Fig. S10. Specifically, in Fig. S10a, equilibrium ternary liquids (at 1500 K) simulated by ab initio MD and EAM MD are compared. In Fig. S10b, a 200-atom Cu46Zr47Al7 MG was first obtained with ab initio MD quenching. The initial liquid structure equilibrated at 1500 K was quenched to 300 K at 1013 K/s, with the density (volume) fixed at 6.926 g/cc [12]. Meanwhile, another 200-atom MG sample was prepared following exactly the same procedure, but using EAM MD. The partial PDFs of the two samples are compared, and the match (within the error and statistical fluctuation) directlydemonstrates the reliability andapplicability of our EAM potentialsin studying the structure of Cu-Zr-Al ternary BMGs. However, for the three reasons discussed at the beginning of the paper, these small samples quenched at an extremely high cooling rate are not suitable for a detailed structural analysis.

Figure S9. Validation of EAM potentials by comparing Cu46Zr54 partial PDFs at different temperatures.(a)Partial PDFs of liquid Cu46Zr54(200 atoms) at 1500 K obtained using EAM (black solid lines) andab initio MD (red dotted lines). The density is fixed at 6.76g/cc, at which the average pressure is below 1 GPa for both EAM and ab initiosimulations.(b) Evolution of the partial PDFs (colored lines) as a function of temperature during NPT quenching (EAM MD, 500 atoms) ofCu46Zr54. The final configuration at 0 K was further relaxed using the conjugated gradient method with ab initio many-body descriptions, and the partial PDFs of thisoptimized structure are also shown(circles).

Figure S10. (a) Partial PDFs of equilibrium liquid Cu46Zr47Al7(200 atoms) at 1500 K obtained using EAM (black solid lines) andab initio MD (red dotted lines).The density is fixed at 6.49g/cc, at which the average pressure is below 1 GPa for both EAM and ab initiosimulations.(b) Comparison of partial PDFs of the ab initio MD-quenched (red dashed lines) and EAM MD-quenched (black solid lines) Cu46Zr47Al7 MGs (200 atoms, quenching rate 1013 K/s). Representative error bars (large statistical fluctuations due to the small sample size and fast quenching rate) for the EAM MD quenching are labeled for Zr-Zr pair, estimated by 10 quenching runs from independent high-temperature liquid structures.

B2.2 Comparisonof glass formation: EAM vs.experimental observation

We also carried out MD simulations to examine the glass forming zone of the Cu-Zr system. To this end, CuxZr100-x alloys (10,000 atoms) with desired compositions were equilibrated at 2000 K. The ensembles (NPT) were then quenched at a constant cooling rate of 1011 K/s, which is comparable with the fastest quenching rate attainable in experiments (e.g. vapor deposition). In Fig. S11, the potential energies as a function of temperature are shown on the left. The final configurations at 300 K were analyzed using common neighbor analysis as well as Voronoi analysis to examine the crystallinity of the samples. It has been found that the glass forming composition at the given cooling condition falls in the range of30 < x <80 (see Fig. S11), which is close to experimental observations [13,14]. This consistency further demonstrates the effectiveness of the EAM potentials developed.Note that the glass forming range may depend on the cooling rate. So this comparison is meant to be only qualitative.