Local structural evolution of Fe₅₄C₁₈Cr₁₆Mo₁₂ bulk metallic glass during tensile deformation and a temperature elevation process: a molecular dynamics study

Abstract

The mechanical and thermal properties of Fe₅₄C₁₈Cr₁₆Mo₁₂ bulk metallic glasses (BMGs) were investigated by a molecular dynamics simulation with the 2NN modified embedded-atom method (MEAM) potential. The fitting process of the cross-element parameters of 2NN MEAM (Fe–C, Fe–Cr, Fe–Mo, C–Cr, C–Mo, and Cr–Mo) was carried out first by the force matching method (FMM) on the basis of the reference data from density functional theory (DFT) calculations. With these fitted parameters, the structure of Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG was constructed by the simulated-annealing basin-hopping (SABH) method, and the angle distribution range of the X-ray diffraction profile of the predicted Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG closely matches that of the experiment profile, indicating the fitted 2NN MEAM parameters can accurately reflect the interatomic interactions of Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG. The Honeycutt–Andersen (HA) index analysis results show a significant percentage of icosahedral-like structures within Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG, which suggests an amorphous state. According to the tensile test results, the estimated Young's modulus of Fe₅₄Cr₁₆Mo₁₂C₁₈ bulk metallic glass is about 139 GPa and the large plastic region of the stress–strain curve shows that the Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG possesses good ductility. Local strain distribution was used to analyze the deformation mechanism, and the results show that a shear band develops homogeneously with the tensile fracture angle (θ_T) at about 50 degrees, in agreement with experimental results 45° < θ_T < 90°. For the temperature elevation results, the discontinuity of the enthalpy–temperature profile indicates the melting point of Fe₅₄Cr₁₆Mo₁₂C₁₈ BMG is about 1310 K. The diffusion coefficients near the melting point were derived by the Einstein equation from the mean-square-displacement (MSD) profiles between 800–1400 K. On the basis of diffusion coefficients at different temperatures, the diffusion barriers of Fe₅₄Cr₁₆Mo₁₂C₁₈ can be determined by the Arrhenius equation. The diffusion barriers of total for Fe, Cr, Mo, C are 31.88, 24.68, 35.26, 22.50 and 31.79 kJ mol⁻¹, respectively. The diffusion barriers of Fe and Cr atoms are relatively lower, indicating Fe and Cr atoms more easily diffuse with the increasing temperature.

---

Introduction

In 1960, Duwez¹ produced metallic glass using fast cooling rates, successfully fabricating an Au–Si system amorphous alloy. The formation of a non-crystalline structure is difficult, however, and in spite of the successful production of this amorphous alloy, there were size limitations. In early 1990, Inoue² proposed an empirical theory regarding the formation of an amorphous alloy that could more easily determine its composition and reduced the importance of cooling rates. As per Inoue's theory, forming an amorphous alloy requires three conditions: (1) a large negative heat, (2) three or more composition elements, and (3) the atomic size between elements must differ by more than 12%. In addition, Inoue also measured the thermal properties of amorphous alloys to investigate their glass-forming ability, and found that they require a crystallization temperature (T_x), glass transition temperature (T_g) and a supercooled liquid region (ΔT_x). With rising interest in amorphous alloys, studies found some elements that possessed high glass-forming ability, such as zirconium, palladium, iron, titanium, and copper. In addition, with advances in technology, the size of amorphous alloys became larger, which meant that bulk metallic glasses (BMG) could be produced. Diameters of bulk metallic glasses have been reported up to 2–30 mm,^3,4 which meant that amorphous alloys could be used in industry or other fields.

Bulk metallic glasses has become a potential material for applications due to their superior mechanical properties, such as high yielding strength up to 6 GPa,^5,6 good ductility⁷ and excellent corrosion resistance.⁸ Accordingly, BMGs have attracted much attention, with some BMGs reported to be suitable for industrial or biomedical applications. For examples, it was reported that Zr-based BMGs present high biocompatibility with tests in NaCl and PBS (phosphate-buffered saline) solutions, showing good corrosion resistance to chloride (Cl) and phosphorus (P) ions.^9–11 Fe-based BMGs display a unique magnetic property, so it is often used in coil production.^12–14 In addition, Fe-based BMGs also have excellent corrosion resistance, allowing for their placement in the body as a biomedical implant material, of which 316L stainless steel, Ti-based BMG, and Co–Cr alloys are commonly used. Although those medical alloys are similar, studies have shown that 316L stainless steel is potentially harmful to the human body because its components contain nickel (Ni), which easily reacts to produce nickel ions after corroding. In recent years, some studies have indicated that the corrosion resistance of Fe-based BMG without nickel is greater than 316L stainless steel.^15,16 In addition, Niinomi's¹⁷ study reported that the Young's modulus of 316L stainless steel (180 GPa) is higher than human bone (about 10–30 GPa), resulting in a stress shielding effect which induces the atrophy of human bone, leading to the loosening of the implanted alloy and refracturing of the human bone.

In 2010, Pilarczyk et al. successfully produced Fe₅₄Cr₁₆Mo₁₂C₁₈ BMG¹⁸ and measured thermal properties such as glass transition temperature and crystallization temperature to investigate the glass forming ability of this BMG at different sizes. From these experimental results, Fe₅₄Cr₁₆Mo₁₂C₁₈ BMG shows a supercooled liquid temperature region arranging between 24 and 79 K, and this phenomena implies Fe₅₄Cr₁₆Mo₁₂C₁₈ BMG possesses a rather high glass forming ability and high thermal stability of the supercooled liquid. However, this study did not provide mechanical properties of Fe₅₄Cr₁₆Mo₁₂C₁₈, and no relevant studies have further reported them. Molecular dynamics simulation of mechanical properties of bulk metallic glass has shown a good degree of accuracy.^19,20 In addition, it is very difficult to investigate the detailed local atomic arrangements around each compositional element and the variations of local atomic arrangements under external loading by the experimental approach directly. The possible alternative to investigate the local atomic arrangement of BMGs is by using numerical simulation. Among different numerical methods, molecular dynamics (MD) simulation can overcome the limitations of traditional empirical approaches and enable detailed observations on local structural variations and the deformation mechanism of BMGs under external loading on the atomic scale.

Therefore, this study will investigate mechanical and structural properties of Fe₅₄Cr₁₆Mo₁₂C₁₈ by molecular dynamics simulation. To the best of our knowledge, this study is the first to provide the interaction parameters between C and three other metal elements for this multi-element system by the force-matching method (FMM).²¹ By these potential parameters, the BMG model was constructed and the detailed local structural arrangements around each atom type were conducted. We also analyzed the changes in atomic structure under the tension test by HA pair analysis,²² and provided some explanations of the simulation results. In the future, we expect that this study can detail the design criterion of new materials for use in biomedical applications.

Simulation model

In order to model the Fe–Cr–Mo–C alloy system by molecular dynamics (MD) simulation, the potential functions, employing the second-nearest neighbor modified embedded-atom method (2NN MEAM),^23,24 were used to describe the interactions between different atomic pairs.

The 2NN MEAM potential form is shown as eqn (1):

where F is the embedding energy which is a function of the atomic electron density , and ϕ_ij is the pair potential interaction. The force-matching method (FMM)²⁵ was used to determine the 2NN MEAM parameters for Fe–C, Fe–Cr, Fe–Mo, Cr–C, Cr–Mo and Mo–C pairs. FMM is based on the variable optimization process of an objective function, and is constructed by the summation of squares of differences between the atomic forces obtained by a potential function and the corresponding atomic forces by ab initio or density functional theory (DFT) calculations. These parameters can be seen in the ESI.†

After the parameters are fitted, they are used to generate the stable Fe₅₄Cr₁₆Mo₁₂C₁₈ amorphous structure by the simulation annealing basin-hopping (SABH) method²⁶ along the search direction for the energy local-minimal structure at higher energy. The unit cell with a total of 4000 atoms (2160 Fe, 640 Cr, 480 Mo and 720 C atoms) is shown in Fig. 1(a), and the model shown in Fig. 1(b) for the tension test by MD was constructed by replicating the unit cell to 6 × 3 × 6 for the x, y and z axes.

Fig. 1 Structures of (a) unit cell and (b) tensile test model for Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG.

Next, the MD simulation was performed by the large-scale atomic/molecular massively parallel simulator (LAMMPS) developed by Plimpton and co-workers.²⁷ By MD simulation, the model was quenched from 1000 K to 300 K for 10 ps to relax the system with an NPT ensemble at 0 GPa. During the tensile process, the periodic boundary conditions (PBC) were applied to the x-, and y-dimensions and the open boundary was used in the z-dimension. The strain rates of 5 × 10⁸ m s⁻¹ were examined to obtain the appropriate strain rate for the current system. During the tension process, the tensile stress at different strains was calculated by the following equation in LAMMPS code:²⁸

where is the effect of the pairwise potential function and indexes 1 and 2 indicate two atoms interacting with each other. The term is the effect of bond angle of the three atoms involved in the interaction. Indexes of m and n represent the m plane and n-direction. The r is the interatomic distance between two atoms 1 and j, m is the weight of the atom, and V_i is a local volume defined by:

Results and discussion

The Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG unit cell of 4000 atoms obtained by SABH method was characterized by the X-ray diffraction (XRD) module REFLEX^29–31 in Accelrys Materials Studio 7.0.³² In REFLEX, Bragg's law is used to obtain the constructive interference intensity for X-rays scattered by materials, and its formula is listed as follows:where θ is a certain angle of incidence when the cleavage faces of crystals appear to reflect X-ray beams. The term d is the distance between atomic layers in a crystal, and λ is the wavelength of the incident X-ray beam. When n is an integer, the diffraction is constructive with higher intensity. While n is a half integer, the diffraction is destructive and the intensity approaches zero.

The XRD and radial distribution function (RDF) profiles of Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG are shown in Fig. 2(a) and (b). One can see no specific crystalline peak appearing in the 2θ range between 20–100° for the XRD profile, and the range of the XRD peak is located at 40–60° and maximum intensity is around 50°, which is consistent with the previous experimental XRD profile of Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG.¹⁸ This indicates our model is in the amorphous state and the fitted cross-element potential parameters of 2NN MEAM can accurately predict the atomic arrangement of Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG as used in the related experiments. For the RDF profile, the broad splitting second peak between 3 and 5 Å indicates the amorphous configuration of Fe₅₄C₁₈Cr₁₆Mo₁₂, which is consistent with the inference of the short range order by the XRD profile. According to the XRD and RDF profiles shown in Fig. 2(a) and (b), the Fe₅₄C₁₈Cr₁₆Mo₁₂ structure constructed by SABH is amorphous and corresponds to realistic Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG in experiment.

Fig. 2 (a) Simulated XRD pattern and (b) RDF profile for the Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG.

A further study into the local microstructural distribution for Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG was conducted by using the Honeycutt–Anderson (HA) pair analysis. The detailed definition of the HA index can be found elsewhere^33–35 and is not presented here. The HA indexes of 1421 and 1422 represent F.C.C. and H.C.P. crystal structures, and 1431, 1541, and 1551 which occupy the largest fraction in the amorphous or liquid state, are used to search the icosahedral local structures. The 1551 pair is particularly characteristic of the icosahedral ordering; the 1541 and 1431 are indexes for the defect icosahedra and F.C.C. defect local (or distorted icosahedra) structures, respectively. HA indexes 1661 and 1441 are employed to identify the local B.C.C. structure. Finally, the indexes 1321 and 1311 are the packing related to rhombohedral pairs which tend to evolve when the 1551 packing forms, which can be viewed as a side product accompanying icosahedral atomic packing. The schematic diagrams for the HA indexes introduced above are illustrated in Fig. 3(a).

Fig. 3 (a) Schematic diagrams corresponding to several characteristic HA indexes; (b) HA index numbers for Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG.

Fig. 3(b) shows the HA index distribution of Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG, and the fraction of icosahedra-like local structures (1551, 1541, and 1431) are about 47.7%. The fraction of the distorted icosahedral structure (1431) occupies the highest fraction of 21.6% among three icosahedra-like fractions, whereas the fraction of perfect icosahedral local structure (1551) is the lowest, with the occupancy of 11.2%. For other HA indexes, the B.C.C. local structures (1441 and 1661), H.C.P. local structure (1422), rhombohedral local structures (1321 and 1311), and F.C.C. local structure (1421) are about 12.1%, 6.9%, 25.3%, and 8.1%, respectively. The high HA fractions of icosahedra-like structures verify the amorphous Fe₅₄C₁₈Cr₁₆Mo₁₂ structure and are consistent with the HA analysis results reported previously for BMGs. Because the second index of the HA analysis is the number of common neighbor atoms between the investigated atomic pair, the higher fractions of 1321 and 1311 HA indexes indicate the local structures of a BMG are more open-packed.³⁶ Compared with the HA analysis of BMGs shown in previous studies,^35,37 the Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG possesses relatively lower icosahedral-like HA fractions and relatively higher rhombohedral HA fractions, indicating more loose local structures distribute within the Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG.

Since the atomic radii of Fe, C, Cr, and Mo are 1.40, 0.70, 1.40, and 1.45 Å, with the atomic size of C smaller than the other three by about 50.0–51.7%, the HA fraction distributions for different atom type pairs are very different. Because the HA index profiles shown in Fig. 3(b) do not contain enough information about the HA fraction distributions for different atom pairs, they should be analyzed to better understand the local structural arrangements around different atom pairs with different pair lengths. Therefore, a more detailed analysis of the HA indexes of different atomic pairs in Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG are shown in Fig. 4. Because Fe occupies the highest atomic fraction in Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG, the Fe-related HA indexes (Fe–Fe, Fe–C, Fe–Cr, and Fe–Mo) are relatively higher than those for other atom type pairs, which can be seen in the first four fractions of each HA index of Fig. 4. The summations of the icosahedra-like HA indexes 1551, 1541, and 1431, referring to the liquid local structures, are about 23.02%, 13.14%, 7.81% and 8.06% for the atom pairs of Fe–Fe, Fe–C, Fe–Cr, and Fe–Mo, respectively. Although the element fraction of Fe is the highest in Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG, one can see the icosahedra-like HA fractions of Fe–Fe and Fe–C are comparable. For the rhombohedral HA fractions, the Fe–Fe pair forms the highest fractions of 1321 and 1311 among all pair types.

Fig. 4 The HA indexes for different pairs of the Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG.

Table 1 lists the average coordination numbers (CNs) of Fe, C, Cr, and Mo atoms in Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG as well as the partial coordination numbers of different atomic pairs. The coordinate number was calculated by counting the amount of first neighbor atoms around the center atom. The cutoff length for the CN calculation was estimated from the first minimal distance of the RDF profile as shown in Fig. 2(b). The first subscript for atomic pair indicates the type of the reference atom and the second subscript stands for the atom type of the first neighbor of the reference atom. Among the average CNs of these four elements, the Cr and Fe atoms have the highest and lowest CNs of 12.48 and 11.70, respectively. A closer investigation of the partial CNs of Fe–Fe, Cr–Fe, Mo–Fe, and C–Fe shows that the partial CN of Cr–Fe is the highest and that of Fe–Fe is the lowest.

Table 1 Average coordination numbers (CNs) for Fe, C, Cr and Mo atoms and each pair in Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG

Fe54C18Cr16Mo12 BMG
Type	Fe–Fe	Fe–Cr	Fe–Mo	Fe–C	Fe_Total
Nij	5.77	2.16	1.48	2.29	11.70
Type	Cr–Fe	Cr–Cr	Cr–Mo	Cr–C	Cr_Total
Nij	6.54	1.97	1.69	2.28	12.48
Type	Mo–Fe	Mo–Cr	Mo–Mo	Mo–C	Mo_Total
Nij	5.99	2.32	1.45	2.40	12.16
Type	C–Fe	C–Cr	C–Mo	C–C	C_Total
Nij	6.20	2.05	1.56	2.13	11.94

---

The Warren–Cowley chemical short-range-order (CSRO) analysis³⁸ for Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG was employed to quantify the attraction and repulsion between element pairs. With the CN information shown in Table 1, the chemical affinities of a referenced atom with its first neighbor atoms are evaluated by the CSRO parameter. The definition of this parameter is as the following equation:

where α_ij is the CSRO parameter of the i-type referenced atom relative to j-type atom, N_ij is the partial CN for the i-type referenced atom relative to j-type atom, and c_j and N_i are the fractions of j-type atom within the alloy and the average CN of i-type atoms, respectively. The value of c_j by N_i is an ideal partial CN for the referenced i-type atom relative to the first neighbor j-type atom, and this value completely depends on the respective atomic composition fraction of Fe₅₄C₁₈Cr₁₆Mo₁₂. According to this analysis, it can be seen that if the alloys are in the ideal solution, the value of N_ij should be very close to c_jN_i and allows the value of α_ij to be 0. On the other hand, a value larger than 0 means the j-type atoms are less prone to gather around i-type atoms, causing the value of N_ij to be less than the number of ideal N_ij. For the same reason, a negative value means the j-type atoms are prone to gather around i-type atoms.

The CSRO parameters of all pairs of Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG are listed in Table 2. The results show that the CSRO parameters for Fe–Fe, Cr–Cr, and Mo–Mo are positive and that of C–C is relatively close to zero. This CSRO analysis result indicates the C atom has no preference to another C atom, and the three other elements display less affinity to themselves, indicating that this alloy easily forms the glass-like structure. Furthermore, most CSRO parameters of C-related pairs are negative except for C–Fe, indicating that the affinities between C and the three other metal elements are relatively higher than Fe-related, Cr-related, and Mo-related ones, which reveals that the smallest atom, C, tends to pair with a metal atom instead of itself.

Table 2 CSRO parameters (α_ij) for all atomic pairs of Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG

Fe54C18Cr16Mo12 BMG
Type	Fe–Fe	Fe–Cr	Fe–Mo	Fe–C
αij	0.086	−0.155	−0.54	−0.085
Type	Cr–Fe	Cr–Cr	Cr–Mo	Cr–C
αij	0.030	0.013	−0.129	−0.016
Type	Mo–Fe	Mo–Cr	Mo–Mo	Mo–C
αij	0.111	−0.163	0.032	−0.067
Type	C–Fe	C–Cr	C–Mo	C–C
αij	0.039	−0.075	−0.088	0.007

---

Fig. 5 shows the stress–strain profile and ΔV/V value with strain for the Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG under tension. One can see the stress increases linearly with strain while strain increases from 0 to 0.05, indicating the elastic behavior of Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG is located within this strain range. The Young's modulus derived from the slope of the stress–strain profile between strains of 0 and 0.02 is about 139 GPa. At strains from 0.05 to 0.1, the stress displays a parabolic increase with increasing strain, and reaches its maximum value at about 6.83 GPa. At strains from 0.1 to 0.5, the stress shows the gradual decrease from its maximal value, indicating the occurrence of fracture. The stress–strain profile for Fe₅₄C₁₈Cr₁₆Mo₁₂ also shows a large plastic region, which can exceed 40%. This result is consistent with those shown in the recent related Fe-based studies.^39,40

Fig. 5 The stress–strain curve and ΔV/V for Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG.

The ΔV/V value, the ratio of open volume (ΔV) to the system volume at strain 0 (V), was used to indicate the volume variation during the tension process. It is a generally held view that more open volume allows for more plastic deformation.⁴¹ The system volume is defined as the summation of the atomic volume calculated by eqn (6), and the ΔV value is calculated by the following equation:

ΔV = V_ε − V₀ (6)

where V_ε and V₀ are the system volumes at strains of ε and 0, respectively. Consequently, the ΔV/V value indicates the percentage increase in system volume during the tension. In Fig. 5, it is apparent that ΔV/V is linearly proportional to the strain within the elastic region from strain 0 to 0.05. At strains from 0.05 to 0.1, the open volume increases parabolically with the strain, and the linear increase of ΔV/V with the strain becomes more significant when the strain exceeds 0.1.

The atomic local shear strain η^Mises_i of an individual atom, introduced by Shimizu et al.,⁴² was used to monitor the development of shear transition zones (STZ) and the formation of the shear band within Fe₅₄C₁₈Cr₁₆Mo₁₂. The detailed definition of η^Mises_i can be found in ref. 43 of this study and is therefore not introduced here. A large η^Mises_i value indicates atom i is under local plastic and shear deformation, whereas a small η^Mises_i value implies atom i undergoes a small amount of movement relative to all its first neighbor atoms or atom i is under local elastic deformation.

Fig. 6(a)–(f) shows the snapshots of Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG with atomic η^Mises_i values at strains of 0, 0.1, 0.2, 0.3, 0.4 and 0.5, which are labelled as letters (a)–(f) on the stress–strain curve of Fig. 5. For the reference structure at strain of 0, the η^Mises_i value of each atom is 0, and the initialization of STZs labeled with black dashed circles in Fig. 6(b) occurs at strain of 0.1. These STZs distribute randomly within Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG. From Fig. 5, one can infer that the sufficient open volume increase significantly activates the initialization of shear banding and enhances the appearance of STZs when the strain exceeds the yielding strain. At strain of 0.2, the extension of STZs begins to form several shear bands, as indicted by the black dashed lines shown in Fig. 6(c), and more shear bands can be seen in Fig. 6(d) at strains of 0.3. In Fig. 5, the ΔV/V value increases more significantly with the strain when strain exceeds 0.1. From the η^Mises_i distributions shown in Fig. 6(c) and (d), one can note that the increase in the shear band number results in the significant increase of the open volume and a local structure rearrangement. The shear bands propagate at a direction 50° from the tensile direction and intersect with one another, resulting in the vein-like pattern. This vein-like pattern can be also seen in previous theoretical⁴⁴ and experimental studies.^45,46 These results show that good ductility of Fe₅₄C₁₈Cr₁₆Mo₁₂ might be caused by the homogeneous development of shear bands which increase the deformation area. Fig. 6(e) and (f) shows the fracture areas at strains of 0.4 and 0.5, where considerable atomic rearrangements occur.

Fig. 6 The distribution of atomic local shear strain.

To understand the local structural rearrangement during the tension process, the numbers of different HA pairs for Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG at different strains during the tension are illustrated in Fig. 7. The vertical axis represents the total number of one particular HA pair and the horizontal axis is the tensile strain. The reason that this study uses the number of pairs instead of pair fraction is that total bonding pairs at each strain will decrease due to increasing distance between atoms during the tensile process such that the fraction is not the best choice to represent the variation of local structure. It can be seen from Fig. 7 that the rhombohedral HA local structure 1311 is predominant at strain higher than 0.05. Among all HA pairs, there are three indexes with notable changes within the strain range 0 to 0.2. The 1551 and 1541 HA indexes, referring to the icosahedral and defected icosahedral structures, show considerable decreases such that the pair numbers of 244 523 and 351 527 for 1551 and 1541 at the beginning of the tensile test decrease to 209 597 and 314 310 at strain of 0.5. For 1311, this HA index increases with the strain from 553 386 to 628 940. Since the second index of the HA analysis is the number of common neighbor atoms between the investigated atomic pair, the decrease of HA indexes with the larger second indexes implies the local structures in Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG become less dense during the tensile process. Note that because the third digit of the HA index represents the bonded number between the first neighbor atoms around the root HA pair, a lower third digit in the HA index indicates a less dense local structure if the first two HA index digits are the same. For example, 1541 is more loose than 1551. Consequently, the 1311 HA index indicates a relatively loose local structure when compared to that of the 1321 HA index. In Fig. 7, the reduced 1551 and 1541 local structures mostly transform into 1311 structures, which is the most loose local structure among those indicated by the HA indexes.

Fig. 7 The HA indexes at different strains for Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG.

The thermal behaviors of Fe₅₄Cr₁₆Mo₁₂C₁₈ were investigated by the MD temperature elevation process, starting from an initial temperature of 300 to 2000 K. During this process, the TtN method,⁴⁷ combining the Nose–Hoover thermostat with the Parrinello–Rahman variable shape size ensemble, was applied to the model shown in Fig. 1(a) with periodic boundary conditions in the x, y, and z dimensions. The TtN method was adopted to maintain a constant temperature and a constant stress of 0 during the temperature elevation simulation. The heating process was conducted by applying a temperature increment of 10 K to the system, which follows a relaxation process of 10 ps before applying the next temperature increment.

Fig. 8 shows the enthalpy profile at different temperatures for Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG during the heating process. The enthalpy value was calculated by taking the average of the enthalpy values of the immediately preceding 3 ps during the relaxation process. The enthalpy is determined from the following equation:

H = U + pV (7)

where U is system energy; p is pressure, and V is the system volume.

Fig. 8 Average enthalpy versus temperature for Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG.

In Fig. 8, the enthalpy is linearly proportional to the increasing temperature from 300 to 1310 K and the discontinuity of this profile appears at 1310 K. When the system temperature is higher than 1310 K, the enthalpy is also linearly proportional to the increasing temperature. Consequently, the melting point (T_m) of Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG is around 1310 K, at which the local structures begin to significantly change, leading to the discontinuity of the enthalpy profile at 1310 K.

The mean-square displacement (MSD) profiles at temperatures ranging from 800 to 1400 K for Fe₅₄C₁₈Cr₁₆Mo₁₂ were used to investigate their dynamical properties. The MSD is defined by a function of time as shown in eqn (8):

where r_i(t) represents the position of atom i at delay time t, and r_i(t₀) indicates the reference position of the corresponding atom at reference time t₀; N represents the total atom number of the investigated system. From Fig. 9, it is clear that the slopes of the MSD profile are generally larger with the increasing temperature.

Fig. 9 The mean-square displacement plots (MSD) of Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG.

It is known that the MSD profile is linear to the delay time over the long-time limit, and thus the diffusion coefficients of Fe₅₄C₁₈Cr₁₆Mo₁₂ can be derived from the slopes of MSD profiles after a longer delay time by the Einstein equation:⁴⁸

where D is the self-diffusion coefficient and N is the number of atoms. The MSD profiles of different elements at different temperatures derived from the Einstein equation for the Fe, C, Cr, and Mo diffusion coefficients of Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG at different temperatures.

In Roy's study, they used the Arrhenius equation to derive the diffusion barrier of Zr and Si atoms on the basis of the calculated diffusion coefficient near the melting points of crystal ZrSi₂.⁴⁹ The formula of the Arrhenius equation for describing the diffusion coefficient at different temperatures is:

where Q is the activation energy, T is the temperature, D₀ is the pre-exponential factor, and R is the Boltzmann constant. To calculate the diffusion barrier, the profiles of ln(D) versus 1/T for total, Fe, Cr, Mo and C atoms are shown in Fig. 10. It can be inferred that the diffusion coefficients of Fe₅₄C₁₈Cr₁₆Mo₁₂ significantly increase with increasing temperature. Note that the ln(D) profile is proportional to the inverse of the temperature, and the diffusion barriers of total, Fe, Cr, Mo and C atoms can be derived by the slopes of the ln(D) profiles. The D₀ values and the diffusion barriers are listed in Table 3. The diffusion barriers of total, Fe, Cr, Mo, C are 31.88, 24.68, 35.26, 22.50 and 31.79 kJ mol⁻¹, respectively. For Fe and Cr atoms, the diffusion barriers are relatively lower, indicating that Fe and Cr atoms more easily diffuse with increasing temperature.

Fig. 10 Diffusion coefficient for Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG: (a) total, (b) Fe, (c) Cr, (d) Mo, and (e) C atoms.

Table 3 The estimated pre-exponential factor (D₀) and activation energy (Q) of Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG

Temperature interval	Type	D0 (m^2 s^−1)	Q (kJ mol^−1)
800–1400 K	Total	2.72 × 10^−10	31.88
	Fe	1.13 × 10^−10	24.68
	C	4.28 × 10^−10	35.26
	Cr	7.70 × 10^−11	22.50
	Mo	2.03 × 10^−10	31.79

---

Conclusion

Molecular dynamics simulations have been conducted to investigate the mechanical properties and thermal properties of Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG by tension and temperature elevation processes, respectively. The unit cell of Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG built by the SABH method is examined by XRD analysis and the XRD profile closely matches that of experimental Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG, indicating Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG from our SABH prediction is totally amorphous and the atomic arrangement of our model is very similar to those used in the experiment. From the HA index analyses, the icosahedra-like structures that indicate liquid local structures appear to be the major local structures in the Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG. The considerable icosahedra-like structures improve the glass forming ability by forming denser and more stable local clusters. Furthermore, chemical affinities are investigated by CSRO parameters. The C-related pairs are nearly all negative except for C–Fe and are relatively higher than Fe-related, Cr-related, and Mo-related pairs, which indicate that the C atom tends to pair with a metal atom instead of itself. This demonstrates that a high diversity of atom size in the solution will advance the formation of the amorphous state.

Based on the stress–strain profile obtained from the tensile test, the predicted Young's modulus is about 139 GPa, which is much less than conventional biomedical implants, such as 316L stainless steel and Co–Cr alloy. In addition, the HA pair analysis of variation in the open volume is also employed to monitor the development of STZ⁵⁰ and the evolution of the shear band. The distributions of the stress–strain curve and open volume with strain show linear increases of stress–strain and open volume–strain curves, which suggest an elastic region. Moreover, stress and open volume increase significantly at STZ initialization stages when the strain exceeds 0.1. This can be attributed to an increase in the number of shear bands, resulting in a significant increase of open volume and the activation of local structural rearrangement. In addition, shear bands measured by the local atomic strain develop along a direction 50° from the tensile direction and indicate good ductility of the Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG.

The temperature, at which the discontinuity of Fe₅₄C₁₈Cr₁₆Mo₁₂ BMG enthalpy–temperature profile during temperature elevation appears, is used to indicate the melting point of about 1310 K. The self-diffusion coefficients of Fe₅₄C₁₈Cr₁₆Mo₁₂ at temperatures near the melting point were calculated by the Einstein equation on the basis of the slopes of the MSD profiles at the long-time limit. On the basis of diffusion coefficients at different temperatures, the diffusion barriers of Fe₅₄Cr₁₆Mo₁₂C₁₈ can be determined by the Arrhenius equation. The diffusion barriers of total, Fe, Cr, Mo, C are 31.88, 24.68, 35.26, 22.50 and 31.79 kJ mol⁻¹, respectively. The diffusion barriers of Fe and Cr atoms are relatively lower, indicating Fe and Cr atoms more easily diffuse with the increasing temperature.

References

1. W. Klement, R. H. Willens and P. Duwez, Non-Crystalline Structure in Solidified Gold–Silicon Alloys, Nature, 1960, 187, 869–870.
2. A. Inoue, High-Strength Bulk Amorphous Alloys with Low Critical Cooling Rates, Mater. Trans., JIM, 1995, 36, 866–875.
3. J. Q. Wang, J. Y. Qin, X. N. Gu, Y. F. Zheng and H. Y. Bai, Bulk Metallic Glasses Based on Ytterbium and Calcium, J. Non-Cryst. Solids, 2011, 3357, 1232–1234.
4. Q. K. Jiang, G. Q. Zhang, L. Yang, X. D. Wang, K. Saksl, H. Franz, R. Wunderlich, H. Fecht and J. Z. Jiang, La-based Bulk Metallic Glasses with Critical Diameter up to 30 mm, Acta Mater., 2007, 55, 4409–4418.
5. J. Wang, R. Li, N. Hua and T. Zhang, Co-based Ternary Bulk Metallic Glasses with Ultrahigh Strength and Plasticity, J. Mater. Res., 2011, 26, 2072.
6. Z. Q. Liu and Z. F. Zhang, Strengthening and Toughening Metallic Glasses: the Elastic Perspectives and Opportunities, J. Appl. Phys., 2014, 115, 163505.
7. W. L. Johnson, Bulk Glass-Forming Metallic Alloys: Science and Technology, MRS Bull., 1999, 24, 42–56.
8. A. Inoue, Stabilization of Metallic Supercooled Liquid and Bulk Amorphous Alloys, Acta Mater., 2000, 48, 279–306.
9. J. J. Oak, D. V. Louzguine-Luzgin and A. Inoue, Fabrication of Ni-free Ti-based Bulk Metallic Glassy Alloy Having Potential for Application as Biomaterial, and Investigation of its Mechanical Properties, Corrosion, and Crystallization Behavior, J. Mater. Res., 2007, 22, 1346–1353.
10. H. F. Li, Y. F. Zheng, F. Xu and J. Z. Jiang, In Vitro Investigation of Novel Ni Free Zr-based Bulk Metallic Glasses as Potential Biomaterials, Mater. Lett., 2012, 75, 74–76.
11. X. N. Gu, N. Li, Y. F. Zheng and L. Q. Ruan, In Vitro Degradation Performance and Biological Response of a Mg–Zn–Zr Alloy, Mater. Sci. Eng., B, 2011, 176, 1778–1784.
12. M. Zhang, A. Wang and B. Shen, Enhancement of Glass-forming Ability of Fe-based Bulk Metallic Glasses with High Saturation Magnetic Flux Density, AIP Adv., 2012, 2, 022169.
13. R. Hasegawa, Applications of Amorphous Magnetic Alloys in Electronic Devices, J. Non-Cryst. Solids, 2001, 287, 405–412.
14. Y. B. Wang, H. F. Li, Y. Cheng, S. C. Wei and Y. F. Zheng, Corrosion Performances of a Nickel-free Fe-based Bulk Metallic Glass in Simulated Body Fluids, Electrochem. Commun., 2009, 11, 2187–2190.
15. H. Zohdi, H. R. Shahverdi and S. M. M. Hadavi, Effect of Nb Addition on Corrosion Behavior of Fe-based Metallic Glasses in Ringer's Solution for Biomedical Applications, Electrochem. Commun., 2011, 13, 840–843.
16. C. Zhang, K. C. Chan, Y. Wu and L. Liu, Pitting Initiation in Fe-based Amorphous Coatings, Acta Mater., 2012, 60, 4152–4159.
17. M. Niinomi and M. Nakai, Titanium-Based Biomaterials for Preventing Stress Shielding Between Implant Devices and Bone, Int. J. Biomater., 2011, 1–10, 2011.
18. W. Pilarczyk, R. Nowosielski and A. Januszka, Structure and Properties of Fe–Cr–Mo–C Bulk Metallic Glasses Obtained by Die Casting Method, Journal of Achievements of Materials and Manufacturing Engineering, 2010, 42, 81–87.
19. X. J. Han and H. Teichler, Liquid-to-glass Transition in Bulk Glass-forming Cu₆₀Ti₂₀Zr₂₀ Alloy by Molecular Dynamics Simulations, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 75, 061501.
20. S. W. Kao, C. C. Huang and T. S. Chin, Simulation of Reduced Glass Transition Temperature of Cu–Zr Alloys by Molecular Dynamics, J. Appl. Phys., 2009, 105, 064913.
21. G. Grochola, S. P. Russo and I. K. Snook, On Fitting a Gold Embedded Atom Method Potential Using the Force Matching Method, J. Chem. Phys., 2005, 123, 204719.
22. J. D. Honeycutt and H. C. Andersen, Molecular Dynamics Study of Melting and Freezing of Small Lennard-Jones Clusters, J. Phys. Chem., 1987, 91, 4950–4963.
23. B.-J. Lee and M. I. Baskes, Second Nearest-Neighbor Modified Embedded-Atom Method Potential, Phys. Rev. B: Condens. Matter Mater. Phys., 2000, 62, 8564–8567.
24. B.-J. Lee, M. I. Baskes, H. Kim and Y. K. Cho, Second Nearest-Neighbor Modified Embedded Atom Method Potentials for bcc Transition Metals, Phys. Rev. B: Condens. Matter Mater. Phys., 2001, 64, 184102.
25. F. Ercolessi and J. B. Adams, Interatomic Potentials from 1st-Principles Calculations – the Force-Matching Method, Europhys. Lett., 1994, 26, 583–588.
26. D. J. Wales and J. P. K. Doye, Global Optimization by Basin-hopping and the Lowest Energy Structures of Lennard-Jones Clusters Containing up to 110 Atoms, J. Phys. Chem. A, 1997, 101, 5111–5116.
27. S. Plimpton, Fast Parallel Algorithms for Short-range Molecular Dynamics, J. Comput. Phys., 1995, 117, 1–19.
28. N. Miyazaki and Y. Shiozaki, Calculation of Mechanical Properties of Solids Using Molecular Dynamics Method, JSME Int. J., Ser. A, 1996, 39, 606–612.
29. G. S. Pawley, Unit-Cell Refinement from Powder Diffraction Scans, J. Appl. Crystallogr., 1981, 14, 357–361.
30. H. M. Rietveld, A Profile Refinement Method for Nuclear and Magnetic Structures, J. Appl. Crystallogr., 1969, 2, 65–71.
31. J. E. Post and D. L. Bish, Rietveld Refinement of Crystal-Structures Using Powder X-Ray Diffraction Data, Rev. Mineral., 1989, 20, 277–308.
32. Materials Studio release 7.0, Accelrys Software Inc., San Diego, USA, 2013.
33. Y. C. Lo, J. C. Huang, S. P. Ju and X. H. Du, Atomic structure evolution of Zr–Ni during severe deformation by HA pair analysis, Phys. Rev. B: Condens. Matter Mater. Phys., 2007, 76, 024103.
34. V. Wessels, A. K. Gangopadhyay, K. K. Sahu, R. W. Hyers, S. M. Canepari, J. R. Rogers, M. J. Kramer, A. I. Goldman, D. Robinson, J. W. Lee, J. R. Morris and K. F. Kelton, Rapid Chemical and Topological Ordering in Supercooled Liquid Cu₄₆Zr₅₄, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 83, 094116.
35. H.-L. Chen, S.-P. Ju, T.-Y. Wu, S.-H. Liu and H.-T. Chen, Investigation of the Mechanical Properties and Local Structural Evolution of Ti₆₀Zr₁₀Ta₁₅Si₁₅ Bulk Metallic Glass during Tensile Deformation: a Molecular Dynamics Study, RSC Adv., 2015, 5, 55383.
36. S.-P. Ju, T.-Y. Wu and S.-H. Liu, Mechanical and Dynamical Behaviors of ZrSi and ZrSi₂ Bulk Metallic Glasses: a Molecular Dynamics Study, J. Appl. Phys., 2015, 117, 105103.
37. S.-P. Ju, H.-H. Huang and T.-Y. Wu, Investigation of the Local Structural Rearrangement of Mg₆₇Zn₂₈Ca₅ Bulk Metallic Glasses During Tensile Deformation: a Molecular Dynamics Study, Comput. Mater. Sci., 2015, 96, 56–62.
38. R. N. Singh and F. Sommer, Segregation and Immiscibility in Liquid Binary Alloys, Rep. Prog. Phys., 1997, 60, 57–150.
39. S. F. Guo, L. Liu, N. Li and Y. Li, Fe-based Bulk Metallic Glass Matrix Composite with Large Plasticity, Scr. Mater., 2010, 62, 329–332.
40. S. F. Guo, J. L. Qiu, P. Yu, S. H. Xie and W. Chen, Fe-based Bulk Metallic Glasses: Brittle or Ductile?, Appl. Phys. Lett., 2014, 105, 161901.
41. X. D. Wang, H. B. Lou, J. Bednarcik, H. Franz, H. W. Sheng, Q. P. Cao and J. Z. Jiang, Structural Evolution in Bulk Metallic Glass under High-Temperature Tension, Appl. Phys. Lett., 2013, 102, 051909.
42. F. Shimizu, S. Ogata and J. Li, Theory of shear banding in metallic glasses and molecular dynamics calculations, Mater. Trans., 2007, 48, 2923–2927.
43. J. T. Wang, P. D. Hodgson, J. D. Zhang, W. Y. Yan and C. H. Yang, Effects of Pores on Shear Bands in Metallic Glasses: a Molecular Dynamics Study, Comput. Mater. Sci., 2010, 50, 211–217.
44. K. Albe, Y. Ritter and D. Şopu, Enhancing the Plasticity of Metallic Glasses: Shear Band Formation, Nanocomposites and Nanoglasses Investigated by Molecular Dynamics Simulations, Mech. Mater., 2013, 67, 94–103.
45. J.-J. Oak and A. Inoue, Attempt to Develop Ti-based Amorphous Alloys for Biomaterials, Mater. Sci. Eng., A, 2007, 449, 220–224.
46. J. S. C. Jang, S. R. Jian, C. F. Chang, L. J. Chang, Y. C. Huang, T. H. Li, J. C. Huang and C. T. Liu, Thermal and Mechanical Properties of the Zr₅₃Cu₃₀Ni₉Al₈ based Bulk Metallic Glass Microalloyed with Silicon, J. Alloys Compd., 2009, 478, 215–219.
47. Y. Qi, T. Cagın, Y. Kimura and W. A. Goddard III, Molecular-dynamics simulations of glass formation and crystallization in binary liquid metals: Cu–Ag and Cu–Ni, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 3527–3533.
48. M. Meunier, Diffusion Coefficients of Small Gas Molecules in Amorphous cis-1,4-polybutadiene Estimated by Molecular Dynamics Simulations, J. Chem. Phys., 2005, 123, 134906.
49. S. Roy and A. Paul, Growth of Hafnium and Zirconium Silicides by Reactive Diffusion, Mater. Chem. Phys., 2014, 143, 1309–1314.
50. H. L. Peng, M. Z. Li and W. H. Wang, Structural Signature of Plastic Deformation in Metallic Glasses, Phys. Rev. Lett., 2011, 106, 135503.

---

Footnotes

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c5ra18168b

‡ Both authors contributed equally to this work.

---