Crystallization of FCC and BCC Liquid Metals Studied by Molecular Dynamics Simulation

: The atomic structure variations on cooling, vitriﬁcation and crystallization processes in liquid metals face centered cubic (FCC) Cu are simulated in the present work in comparison with body centered cubic (BCC) Fe. The process is done on continuous cooling and isothermal annealing using a classical molecular-dynamics computer simulation procedure with an embedded-atom method potential at constant pressure. The structural changes are monitored with direct structure observation in the simulation cells containing from about 100 k to 1 M atoms. The crystallization process is analyzed under isothermal conditions by monitoring density and energy variation as a function of time. A common-neighbor cluster analysis is performed. The results of thermodynamic calculations on estimating the energy barrier for crystal nucleation and a critical nucleus size are compared with those obtained from simulation. The di ﬀ erences in crystallization of an FCC and a BCC metal are discussed. to


Introduction
Crystallization of metals and alloys from the melt (liquid state) takes place by nucleation and growth [1]. This process has occupied the minds of scientists for quite some time owing to still unresolved issues, especially related to high nucleation and growth rates of pure substances. Pure metals are known to exhibit a high crystal nucleation rate at a low enough temperature [1] and high growth rates in a wide temperature range [2]. Crystallization of different substances [3,4] including liquid metals [5] and alloys [6] has been studied by molecular-dynamics (MD) simulation with various potentials. Furthermore, MD simulation of glass transition of pure metals, such as Fe [7,8], Ni [9,10], Cu [11], Al [12] and other metals [13] has been performed. In a recent work possible heterogeneous nucleation of Fe was studied at 0.67T m and 0.58T m (T m is the melting point which in pure metals is equal to the liquidus temperature T l ) [14]. However, one should mention that T m of 2400 K given by the potential used is significantly overestimated compared to the experimental value. Recent embedded atom potentials provide a better correspondence to the experimental values. The crystallization behavior of other metals was also studied via MD simulation [15] and compared with the experimental data [16]. The energy barriers for crystallization of pure metals were also compared with those obtained from computer simulation [17]. In the present work the atomic structure variation, vitrification and crystallization processes in liquid metals (FCC Cu and BCC Fe) are simulated on continuous cooling and isothermal atomic structure variation, vitrification and crystallization processes in liquid metals (FCC Cu and BCC Fe) are simulated on continuous cooling and isothermal annealing using a classical MD computer simulation procedure with an embedded-atom method potential at constant pressure.

Computational Procedure
Molecular dynamics (MD) computer simulation using a software package for classical molecular dynamics (LAMMPS) [18] was used to model the crystallization process of Fe and Cu metals at periodic boundary conditions. The simulation was performed at 1 fs time step using the embedded atom potentials for Cu [19] and Fe [20] at nearly constant temperature and pressure. In order to let crystallization occur, the cubic cell size was typically chosen to be about 10 nm. A typical crystalline cell of the atoms containing 108,000 for FCC and 128,000 atoms for BCC structure was heated to of 2500 K to melt and then cooled down to a certain temperature. In some cases, larger sized cells were used as indicated. Melting is confirmed by the radial distribution function and stabilization of the density variation with time. Cooling was performed to different temperatures at about 5 × 10 13 K/s. Temperature upon simulation was maintained at ±6 K ( Figure 1) while pressure was kept around zero. Neither Fe nor Cu evaporated in the molten state as should have happened at zero pressure, likely because there was not enough time for a gaseous phase to nucleate from inside the atomic cell (a gaseous phase has to overcome an energy barrier to be formed as well as a crystalline one) and absence of the surfaces owing to periodic boundary conditions. A thermostat was used to control the temperature [21,22] while pressure was maintained by a Barostat [23]. A software package "OVITO" [24] was used to visualize and analyze the simulation results. Adaptive Common Neighbor Analysis (CNA) [25] was used to analyze structural changes in Fe and Cu.

Results and Discussion
The density changes on heating and cooling as a function of temperature are illustrated in Figure 2a. Equilibrium liquid state was obtained at 2500 K in less than 1 ps. The holding time was also varied from 1 to 100 ps without visible changes in the liquid structure and density. Thus, the atomic cells have been kept for at least 10 ps at 2500 K prior to cooling. Neither Fe nor Cu evaporated in the molten state as should have happened at zero pressure, likely because there was not enough time for a gaseous phase to nucleate from inside the atomic cell (a gaseous phase has to overcome an energy barrier to be formed as well as a crystalline one) and absence of the surfaces owing to periodic boundary conditions. A thermostat was used to control the temperature [21,22] while pressure was maintained by a Barostat [23]. A software package "OVITO" [24] was used to visualize and analyze the simulation results. Adaptive Common Neighbor Analysis (CNA) [25] was used to analyze structural changes in Fe and Cu.

Results and Discussion
The density changes on heating and cooling as a function of temperature are illustrated in Figure 2a. Equilibrium liquid state was obtained at 2500 K in less than 1 ps. The holding time was also varied from 1 to 100 ps without visible changes in the liquid structure and density. Thus, the atomic cells have been kept for at least 10 ps at 2500 K prior to cooling. Vitrification of both Fe and Cu liquids was attained at the cooling rate of about 5 × 10 13 K/s (Figure 2). At the cooling rate of 5 × 10 12 K/s and lower the studied metals start to crystallize directly  Vitrification of both Fe and Cu liquids was attained at the cooling rate of about 5 × 10 13 K/s (Figure 2). At the cooling rate of 5 × 10 12 K/s and lower the studied metals start to crystallize directly on cooling. Similarly, high critical cooling rate values indicating high instability of the supercooled liquid were obtained in other works [13,26,27]. However, theoretical calculations [28] and the experimental results [29] suggest a lower critical cooling rate (of about 10 9 K/s) for glass-transition in pure metals, especially with BCC lattice. The difference could be connected with the type of potentials used.
The values of glass-transition temperature (T g ) defined by the change of the density as a function of temperature slope are about 1100 K for Fe and 800 K for Cu (Figure 2a). Of course, these values correspond to this very high cooling rate used when viscosity is still low. However, the value for Fe corresponds quite well to the isochoric conditions when the volume of the liquid becomes similar to that of the crystalline phase [30]. Two typical pair distribution functions (PDF) for Fe are shown in Figure 2b. The values of T g derived from PDF min [31] and ratios of PDF min /PDF max [32] are shown in Figure 2c. These values are about 100 K lower than those obtained from Figure 2a. As one can see for Fe, having a BCC lattice at room temperature has much lower PDF min value compared to Cu which has an FCC lattice. It resembles the crystalline structure: there is a deeper gap between the coordination shells in BCC lattice. Splitting of the first and second pair distribution function (PDF(R)) maxima into two peaks is clearly observed in Figure 2b. It indicates short-and medium-range ordering of the liquid and glassy phases on cooling and that glassy Cu and Fe inherit the short-range order of the corresponding crystals.
Phase transformation kinetics for Fe and Cu are shown in Figure 3; Figure 4, respectively, by monitoring the density changes. Energy of the system changes in the way opposite to density. Please, note that the final product density is different from case to case owing to different fractions of the internal defects. Full system equilibration upon crystallization requires much longer time. While at 1100 K and above the density, volume and energy values are visibly stable within the incubation period (relaxation proceeds in 10 ps for liquid Fe at 1100 K) (Figure 3a), these values gradually change at a lower temperature owing to structural relaxation at low atomic mobility. The structure relaxation time becomes longer than the simulation time below 1100 K for Fe which corresponds quite well to the equivolume/isochoric glass-transition temperature determined before [30]. It also takes place at about 800 K for Cu when the density and energy of the system gradually change with time prior to crystallization. Both BCC and FCC crystals nucleate directly from the melt and no two-stage nucleation [33] is found.
Metals 2020, 10, x FOR PEER REVIEW 4 of 11 on cooling. Similarly, high critical cooling rate values indicating high instability of the supercooled liquid were obtained in other works [13,26,27]. However, theoretical calculations [28] and the experimental results [29] suggest a lower critical cooling rate (of about 10 9 K/s) for glass-transition in pure metals, especially with BCC lattice. The difference could be connected with the type of potentials used. The values of glass-transition temperature (Tg) defined by the change of the density as a function of temperature slope are about 1100 K for Fe and 800 K for Cu ( Figure 2a). Of course, these values correspond to this very high cooling rate used when viscosity is still low. However, the value for Fe corresponds quite well to the isochoric conditions when the volume of the liquid becomes similar to that of the crystalline phase [30]. Two typical pair distribution functions (PDF) for Fe are shown in Figure 2b. The values of Tg derived from PDFmin [31] and ratios of PDFmin/PDFmax [32] are shown in Figure 2c. These values are about 100 K lower than those obtained from Figure 2a. As one can see for Fe, having a BCC lattice at room temperature has much lower PDFmin value compared to Cu which has an FCC lattice. It resembles the crystalline structure: there is a deeper gap between the coordination shells in BCC lattice. Splitting of the first and second pair distribution function (PDF(R)) maxima into two peaks is clearly observed in Figure 2b. It indicates short-and medium-range ordering of the liquid and glassy phases on cooling and that glassy Cu and Fe inherit the short-range order of the corresponding crystals.
Phase transformation kinetics for Fe and Cu are shown in Figure 3; Figure 4, respectively, by monitoring the density changes. Energy of the system changes in the way opposite to density. Please, note that the final product density is different from case to case owing to different fractions of the internal defects. Full system equilibration upon crystallization requires much longer time. While at 1100 K and above the density, volume and energy values are visibly stable within the incubation period (relaxation proceeds in 10 ps for liquid Fe at 1100 K) (Figure 3a), these values gradually change at a lower temperature owing to structural relaxation at low atomic mobility. The structure relaxation time becomes longer than the simulation time below 1100 K for Fe which corresponds quite well to the equivolume/isochoric glass-transition temperature determined before [30]. It also takes place at about 800 K for Cu when the density and energy of the system gradually change with time prior to crystallization. Both BCC and FCC crystals nucleate directly from the melt and no two-stage nucleation [33]      As can be found for pure Fe in Figure 3, the incubation time (defined as an intersection of two tangents to the plot before and after the onset time) changes irregularly with the atomic cell size owing to stochastic character of nucleation. In a cell containing 432,000 atoms, two overcritical nuclei nucleated at about 120 ps. In a cell containing 1,024,000 atoms, the first nucleus was formed at 150 ps, the second at 170 ps, the third at 210 ps, the fourth at 230 ps and the fifth at 240 ps. If one assumes total nucleation time of 90 ps for four nuclei per a cell of 23 nm size, then the average homogeneous nucleation rate is 3.7 × 10 33 m −3 ·s −1 . It is a high value but quite close to that reported in an earlier work [8]. The growing nanocrystals have irregular but nearly spherical shape. In a cell containing 128,000 atoms, the first overcritical stably growing nucleus was formed at 280 ps and consumed the entire volume at about 410 ps. Thus, at these temperatures (1100 and 900 K which are 0.61 Tl and 0.50 Tl, respectively; when the liquidus temperature Tl = 1811 K) crystallization of Fe is hardly statistically predictable, and it is not so sensitive to the cell size if its length exceeds about 10 nm. BCC Cu at 750 K (which is 0.55 Tl) has even higher homogeneous nucleation rate of 8 × 10 34 m −3 s −1 . At the homological temperature higher than 0.65 Tl crystallization of Cu and Fe is not detectable within 3 ns of maximum simulation time.
Owing to rather high nucleation rates, the isothermal crystallization curves can be statistically well reproduced below 0.55 Tl for Cu ( Figure 4b) and below 0.4 Tl for Fe at the chosen atomic cells. An average incubation time is about 2 times longer than that on pure Ni studied earlier [34]. From this viewpoint Ni is a special liquid which is the least stable against crystallization among these metals. Furthermore, in spite of the fact that the nucleation barrier is quite high in pure Ni the nucleation rate is also high attributed to the high atom attachment kinetics [35].
From the isothermal crystallization plots, the beginning, also called an incubation time, and finish time of transformation can be calculated using two corresponding tangents to the plot before and after the inflection point. Note that actual transformation time is longer because the two tails of the S-curve are not taken into consideration when two tangents are applied. However, it becomes more difficult to detect the incubation period for each plot below 700 K owing to the long relaxation time manifested in constant variation of the baseline. Nevertheless, a time-temperaturetransformation (TTT) diagram was constructed, and it is shown in Figure 5. The simulation results are well reproduced (owing to high nucleation rate) below 0.55 Tl for Cu and below 0.4 Tl for Fe. The nose of the TTT diagram is about 0.55 Tl for Fe and about 0.6 Tl for Cu. As can be found for pure Fe in Figure 3, the incubation time (defined as an intersection of two tangents to the plot before and after the onset time) changes irregularly with the atomic cell size owing to stochastic character of nucleation. In a cell containing 432,000 atoms, two overcritical nuclei nucleated at about 120 ps. In a cell containing 1,024,000 atoms, the first nucleus was formed at 150 ps, the second at 170 ps, the third at 210 ps, the fourth at 230 ps and the fifth at 240 ps. If one assumes total nucleation time of 90 ps for four nuclei per a cell of 23 nm size, then the average homogeneous nucleation rate is 3.7 × 10 33 m −3 ·s −1 . It is a high value but quite close to that reported in an earlier work [8]. The growing nanocrystals have irregular but nearly spherical shape. In a cell containing 128,000 atoms, the first overcritical stably growing nucleus was formed at 280 ps and consumed the entire volume at about 410 ps. Thus, at these temperatures (1100 and 900 K which are 0.61 T l and 0.50 T l , respectively; when the liquidus temperature T l = 1811 K) crystallization of Fe is hardly statistically predictable, and it is not so sensitive to the cell size if its length exceeds about 10 nm. BCC Cu at 750 K (which is 0.55 T l ) has even higher homogeneous nucleation rate of 8 × 10 34 m −3 s −1 . At the homological temperature higher than 0.65 T l crystallization of Cu and Fe is not detectable within 3 ns of maximum simulation time.
Owing to rather high nucleation rates, the isothermal crystallization curves can be statistically well reproduced below 0.55 T l for Cu ( Figure 4b) and below 0.4 T l for Fe at the chosen atomic cells. An average incubation time is about 2 times longer than that on pure Ni studied earlier [34]. From this viewpoint Ni is a special liquid which is the least stable against crystallization among these metals. Furthermore, in spite of the fact that the nucleation barrier is quite high in pure Ni the nucleation rate is also high attributed to the high atom attachment kinetics [35].
From the isothermal crystallization plots, the beginning, also called an incubation time, and finish time of transformation can be calculated using two corresponding tangents to the plot before and after the inflection point. Note that actual transformation time is longer because the two tails of the S-curve are not taken into consideration when two tangents are applied. However, it becomes more difficult to detect the incubation period for each plot below 700 K owing to the long relaxation time manifested in constant variation of the baseline. Nevertheless, a time-temperature-transformation (TTT) diagram was constructed, and it is shown in Figure 5. The simulation results are well reproduced (owing to high nucleation rate) below 0.55 T l for Cu and below 0.4 T l for Fe. The nose of the TTT diagram is about 0.55 T l for Fe and about 0.6 T l for Cu. According to adaptive common neighbor analysis (Figure 6), crystallization becomes detectable by volume and energy changes when the fraction of atoms in crystalline clusters attains about 2%. The increase in the numbers of atoms with the icosahedral-like configuration prior to the nucleation of crystalline phases found in earlier works [36,37] is not seen in the present study (see Figure 6). Although only BCC crystals form and grow in Fe FCC, BCC and hexagonal close packed (HCP) atomic arrangements are found in Cu owing to low stacking fault energy of this metal (well estimated experimentally) as seen in Figure 7.  According to adaptive common neighbor analysis (Figure 6), crystallization becomes detectable by volume and energy changes when the fraction of atoms in crystalline clusters attains about 2%. The increase in the numbers of atoms with the icosahedral-like configuration prior to the nucleation of crystalline phases found in earlier works [36,37] is not seen in the present study (see Figure 6). Although only BCC crystals form and grow in Fe FCC, BCC and hexagonal close packed (HCP) atomic arrangements are found in Cu owing to low stacking fault energy of this metal (well estimated experimentally) as seen in Figure 7. According to adaptive common neighbor analysis (Figure 6), crystallization becomes detectable by volume and energy changes when the fraction of atoms in crystalline clusters attains about 2%. The increase in the numbers of atoms with the icosahedral-like configuration prior to the nucleation of crystalline phases found in earlier works [36,37] is not seen in the present study (see Figure 6). Although only BCC crystals form and grow in Fe FCC, BCC and hexagonal close packed (HCP) atomic arrangements are found in Cu owing to low stacking fault energy of this metal (well estimated experimentally) as seen in Figure 7.  Typical growth rates (crystal radius change as a function of time) are shown in Table 1. The growth rates are lower than those found earlier in simulations [15]. The growth rate of FCC Cu crystals is less temperature dependent than that of BCC Fe. Furthermore, one should note that some growing particles have much lower growth rate owing to interaction with other nuclei (even undercritical) [38]. At the same time the observed growth rate values for Fe are quite close to the experimental ones of 40 m/s [39]. The observed growth rate of the order of tens meters per second are lower than those obtained for flat interfaces. A liquid/solid interfacial energy (σ) value of 0.37 J/m 2 is obtained in the present work for Fe at 1100 K as an excess energy per surface area of a single particle of 1.8 nm radius. It is calculated taking into account the potential energy difference of 23.7 eV at 170 ps of the simulation time between its value for the entire atomic cell and a sum of that energy values for the atoms belonging to a crystalline particle of about 3.6 nm diameter formed at that stage (1657 atoms) and the potential energy of the other atoms in the cell belonging to liquid phase (128,000−1657 = 126,343 atoms). The potential energy for the liquid part of the atomic cell was obtained at the 50 ps stage when all atoms belonged to the supercooled liquid phase while that of a crystal was taken for fully crystalline state on heating the initial perfect crystal to 1100 K. The value of 0.33 J/m 2 for the Cu crystal at 850 K is found in a similar way. These values are not too far from the experimental values of 0.28 and 0.18 J/m 2 obtained for Fe and Cu, respectively [1]. Please note that although the calculated values are somewhat overestimated owing to the internal defects in the growing crystals, the value for Fe is still higher than that for Cu. The earlier simulation results indicated that the crystal-melt interfacial energy calculated from a planar interface needs to be adjusted to a smaller value in order to predict the correct nucleation rate at deep supercooling [40].
Calculation of the critical nucleus size (Rc) and work for the formation of critical nucleus (W*) was done according to classical nucleation theory: Typical growth rates (crystal radius change as a function of time) are shown in Table 1. The growth rates are lower than those found earlier in simulations [15]. The growth rate of FCC Cu crystals is less temperature dependent than that of BCC Fe. Furthermore, one should note that some growing particles have much lower growth rate owing to interaction with other nuclei (even under-critical) [38]. At the same time the observed growth rate values for Fe are quite close to the experimental ones of 40 m/s [39]. The observed growth rate of the order of tens meters per second are lower than those obtained for flat interfaces. A liquid/solid interfacial energy (σ) value of 0.37 J/m 2 is obtained in the present work for Fe at 1100 K as an excess energy per surface area of a single particle of 1.8 nm radius. It is calculated taking into account the potential energy difference of 23.7 eV at 170 ps of the simulation time between its value for the entire atomic cell and a sum of that energy values for the atoms belonging to a crystalline particle of about 3.6 nm diameter formed at that stage (1657 atoms) and the potential energy of the other atoms in the cell belonging to liquid phase (128,000−1657 = 126,343 atoms). The potential energy for the liquid part of the atomic cell was obtained at the 50 ps stage when all atoms belonged to the supercooled liquid phase while that of a crystal was taken for fully crystalline state on heating the initial perfect crystal to 1100 K. The value of 0.33 J/m 2 for the Cu crystal at 850 K is found in a similar way. These values are not too far from the experimental values of 0.28 and 0.18 J/m 2 obtained for Fe and Cu, respectively [1]. Please note that although the calculated values are somewhat overestimated owing to the internal defects in the growing crystals, the value for Fe is still higher than that for Cu. The earlier simulation results indicated that the crystal-melt interfacial energy calculated from a planar interface needs to be adjusted to a smaller value in order to predict the correct nucleation rate at deep supercooling [40]. Calculation of the critical nucleus size (R c ) and work for the formation of critical nucleus (W*) was done according to classical nucleation theory: Using the Gibbs free energy difference (∆G v ). σ values of 0.277 J/m 2 for Fe and 0.178 J/m 2 for Cu were obtained from Ref. [1]. The temperature dependences of enthalpies of liquid and crystalline phases were used to estimate the enthalpy of melting at about 1809 K. The enthalpy difference of 0.164148 eV/at or 15.837 kJ/mol found for Fe is not far from the experimental value of 15.2 kJ/mol [1]. The value of 13 kJ/mol was used for Cu [1]. ∆G v can be estimated from the liquid crystallization enthalpy (∆H c ) and the liquidus temperatures T l using the following equation [41]: According to thermodynamic calculations the critical nucleus radius in the studied temperature range is about 0.4-0.6 nm for Cu and for 0.5-0.6 nm Fe. Steadily growing particles found in this simulation were larger than 0.8 nm in size. At such small size the energy barrier W* to the thermal energy RT ratio is about 20 for Cu and about 30 for Fe. A higher energy barrier found for Fe is rather responsible for a lower number of nuclei found for Fe at the same homological temperature. It is interesting to note that even at such high energy barriers homogeneous nucleation takes place after hundreds of picoseconds.
The number of crystal nuclei formed within the typical cells of about 10 nm length in all metals ranges from 1 to 10, leading to the number density of precipitates of the order of 10 24 m −3 . The highest number density of nuclei is observed in Cu which leads to better reproducibility of the crystallization traces at different simulation runs. Such high number densities of precipitates are more or less typical for crystallization of pure metals at extremely low homological temperatures of about 0.55 T l or lower, which are unattainable in the experiments on undercooling but quite common in various metallic glasses which nanocrystallize in a similar temperature range [42,43].

Conclusions
Liquid structure variation of cooling and glass transition indicated that glassy Cu and Fe inherit the short-range order of the corresponding crystals. Nevertheless, a nucleation and growth type (with an energy barrier) phase transformation is observed. The crystallization kinetics of FCC Cu and BCC Fe were analyzed in isothermal conditions by monitoring the density and energy variation as a function of time, and a TTT diagram was constructed. BCC Fe shows a lower number density of precipitates than FCC Cu. A higher energy barrier for nucleation found for Fe is rather responsible for a lower number of nuclei found for Fe at the same homological temperature. On the other hand, BCC Fe, in general, shows a shorter incubation period than FCC Cu and generally a higher growth rate though the nucleation rate of Cu crystals at a higher rate. The nose of the TTT diagram is about 0.55 T l for Fe and about 0.6 T l for Cu. Homogeneous nucleation of both metals becomes difficult to detect above about 0.65 T l within a reasonable simulation time scale, while below~0.6 T l a high population density of crystalline precipitates is obtained. The growing nanocrystals have an irregular but nearly spherical shape. The observed growth rates of the order of tens meters per second are lower than those obtained for flat interfaces. Some of the crystalline particles displayed a lower growth rate than the others owing to the lattice defects within the growing crystals and existing order in the supercooled liquid. The crystal growth rate was scattered from crystal to crystal likely owing to the local order in surrounding environment and the lattice defects (mostly HCP staking faults) within the growing crystals in Cu. According to thermodynamic calculations the critical nucleus radius in the studied