Intermetallics 19 (2011) 630e635
Contents lists available at ScienceDirect
Intermetallics journal homepage: www.elsevier.com/locate/intermet
Atomistic structural evolution with cooling rates during the solidification of liquid nickel F. Li a, b, X.J. Liu b, *, H.Y. Hou a, c, G. Chen a, c, G.L. Chen a, b, ** a
Engineering Research Center of Materials Behavior and Design, Ministry of Education, Nanjing University of Science and Technology, Nanjing 210094, China State Key Laboratory for Advanced Metals and Materials, University of Science and Technology Beijing, Beijing 100083, China c Jiangsu Institute of Advanced Materials, Danyang 212300, Jiangsu, China b
a r t i c l e i n f o
a b s t r a c t
Article history: Received 15 October 2010 Received in revised form 24 December 2010 Accepted 28 December 2010 Available online 22 January 2011
In this paper, molecular dynamics simulations based on the second-moment approximation of tightbinding scheme have been performed to investigate the effect of cooling rate on solidification microstructures of liquid Ni. Two transitional structures characterized by one-dimensional (1D) and 2D periodicities respectively between the crystal and amorphous states have been found as the cooling rates range from 6 1012 K/s to 1 1013 K/s. As the cooling rates Q 1.5 1013 K/s, an amorphous structure can be obtained, whilst crystal structures are formed when Q 4 1012 K/s. Moreover, our results reveal that the intensity ratio (g21/g22) of the two subpeaks, which split from the second peak on the pair distribution function for the amorphous state, can act as a structural indicator to differentiate amorphous, translational structure and crystalline states. As such, one may determine the structure state of a material by estimating the value of g21/g22 from its pair distribution function. Ó 2010 Elsevier Ltd. All rights reserved.
Keywords: B. Glasses, metallic D. Microstructure E. Simulations, atomistic
1. Introduction The liquid-to-solid phase transformation is a fundamental issue in the field of condensed matter physics and materials science [1e6]. The solidification phase transformation in which a liquid is transformed into a solid upon cooling generally proceeds in two different ways associated with cooling rates. One route is the crystallization which takes place under a relatively low cooling rate and finally produces crystalline solids with three-dimensional (3D) periodicity. On the other hand, when the cooling rate is large enough to suppress the nucleation and growth of crystalline phases, an amorphization transformation occurs, which leads to form an amorphous state without long-range order. In particular, bulk metallic glasses (BMGs) resulting from the amorphization transformation have attracted considerable attention from both scientific and industrial fields recently, due to their unique properties and potentials for engineering applications [7e11]. Although there are numerous investigations on solidification of metallic liquids upon cooling, the evolution of atomistic structures with cooling
* Corresponding author. ** Corresponding author. State Key Laboratory for Advanced Metals and Materials, University of Science and Technology Beijing, Beijing 100083, China. E-mail addresses:
[email protected] (X.J. Liu),
[email protected] (G.L. Chen). 0966-9795/$ e see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.intermet.2010.12.016
rates during solidification of metallic liquids is far from fully understood. The cooling rate is a very critical factor influencing the final solidification microstructures. In previous experimental studies about rapid quenching of alloys [12e16], they mainly focused on the critical cooling rates for forming glasses to estimate the glassforming ability of these alloys. However, there are few systematic investigations on the evolution of final solidification microstructures with cooling rates, especially the structural evolution at atomic-scale level, due to the limit of experimental techniques. Experimentally, it is still a significant challenge to survey the atomic configuration in an imperfect crystal within 1e2 nm although the advanced high resolution transmission electron microscopy can probe the atomic structure down to sub-angstrom scale. On the other hand, molecular dynamics (MD) simulation is currently one of the most powerful tools for studying the atomic-scale structure as a function of various parameters such as temperature and cooling rate. Furthermore, MD simulations can provide some important information which is beneficial to understanding the experimental observations but difficult to access by experiments themselves. Actually, MD simulations have been successfully used to study the structural, thermodynamic and dynamic properties of the glass formation, crystallization and relaxation process from liquid by using different types of interatomic potentials [17e24]. Among these MD simulations, much attention has been paid to study the process of glass transition or nucleation and growth.
F. Li et al. / Intermetallics 19 (2011) 630e635
Nevertheless, the relationship between the final solidification microstructures and quenching rates, as well as the critical cooling rates for glass formation and crystallization in pure Ni system, is less concerned. In this study, we systematically investigated the relationship between the final solidification microstructures of Ni at atomicscale level and the cooling rate using MD simulations based on the second-moment approximation of tight-binding scheme (TB-SMA). The critical cooling rates for forming glassy and crystalline Ni were determined separately, and the atomic structural evolution from the glass to the crystal with the decrease of cooling rates was also explored. Furthermore, the pair distribution function (PDF) as a function of cooling rate was analyzed, in particular we concerned how the splitting second peak in the PDF of the amorphous state decompose into two separate peaks in the crystal. The result indicates that the intensity ratio of these two peaks can be used to character the structure state. 2. Simulation method TB-SMA is a well-studied semiempirical many-body interaction potential and has been widely used in numerical simulations in metals and alloys. In this study, we adopted the TB-SMA potential fitted by Cleri and Rosato to describe the interatomic interactions between Ni atoms, in which the total potential energy of the system is given by [25]
U ¼
8 N < X N X i¼1
:
j ¼ 1si
N X
Q1 ¼ 5 1013 K/s, Q2 ¼ 2 1013 K/s, Q3 ¼ 1.5 1013 K/s, Q4 ¼ 1 1013 K/ s, Q5 ¼ 8 1012 K/s, Q6 ¼ 6 1012 K/s, Q7 ¼ 4 1012 K/s, Q8 ¼ 1 1012 K/ s, and Q9 ¼ 5 1011 K/s. 3. Results and discussions 3.1. Potential energy and volume Fig. 1 illustrates the average potential energy (Fig. 1a) and molar volume (Fig. 1b) as functions of temperature at various quenching rates. These nine curves can be classified generally into three kinds according to their features. For the cooling rates not less than 1.5 1013 K/s (Q1, Q2 and Q3), only a continuously subtle change in slope of energy and volume with temperature can be seen, showing the characteristics of liquid-to-glass phase transformation. For quenching rates below 6 1012 K/s (Q7, Q8 and Q9), however, an abrupt energy and volume jump with decreasing temperature, which represents the typical liquid-to-crystal first-order phase transformation. As for the intermediate cooling rates (Q4, Q5 and Q6), it is found that the variations of potential energy and molar volume with temperature exhibit an interesting feature which is distinct from those observed in the aforementioned cases. There are no clear inflexions in the temperature-dependent energy and volume curves, but exist a nearly linear relationship between the energy/volume and temperature. This indicates that some transitional structures are formed to bridge the gap between the amorphous and the crystal
rij Aexp p 1 r0
3 9 1=2 = r x2 exp 2q ij 1 5 4 ; r0 j ¼ 1si 2
631
ð1Þ
where rij is the interatomic distance between the atom i and its neighbor j, and r0 is the equilibrium distance between atoms in the perfect crystal. The parameters for Ni used in our simulations are as following [25]: A ¼ 0.0376 eV, x ¼ 1.070 eV, p ¼ 16.999 and q ¼ 1.189. The cutoff radius rc is fixed to be 0.50 nm to ensure the rationality of our simulation results. The simulations were performed using the classical MD program XMD 2.5.32 [26]. In our MD simulations, a system with 4000 particles distributed in a 10 10 10 lattice unit cubic box under periodic boundary conditions was studied. The initial cubic box size was 3.523 nm and the initial positions of atoms were set according to the structure of face-centered cubic (fcc) Ni crystal. The MD time step was chosen as 1 fs to integrate Newton’s equations of motion. The present MD simulations have been carried out as an NPT isobariceisothermal ensemble with constant particle number (N), pressure (P) and temperature (T). This method allows one to study phase transition while the volume of the simulated system is allowed to adjust itself to sustain the given pressure and temperature. The MD simulations were implemented with 5th Gear predictionecorrection algorithm, using the “Pressure Clamp” method to maintain a constant P and the “Temperature Clamp” method to maintain a constant T. Firstly, the well-equilibrated initial system was prepared by gradually heating the perfect crystal from 300 K to 2000 K (the melting point of Ni is 1728 K [27]) at a heating rate of 1.0 1012 K/s under a zero external pressure. And then let the system run 30,000 MD time steps (30 ps) at 2000 K to obtain an equilibrium liquid state, and finally the liquid system was quenched to T ¼ 300 K. During the relaxation process at each given temperature, 100 configurations were saved in order to record volume, energy and position information for analyzing the thermodynamic and structural properties. We simulated the cooling procedure at nine cooling rates, ranging from 1011 to 1013 K/s, i.e.,
Fig. 1. The potential energy (a) and molar volume (b) as functions of temperature for different cooling rates of Ni simulated system.
632
F. Li et al. / Intermetallics 19 (2011) 630e635
states under the modest cooling rates. The instantaneous configuration images corresponding to the transitional structures will be given in the next section and confirm this point. Further analyses show that the potential energy and molar volume at 300 K are respectively 4.06 eV/atom and 7.18 cm3/mol for the amorphous state, 4.09 eV/atom and 7.07 cm3/mol for the transitional state, and 4.14 eV/atom and 6.91 cm3/mol for crystalline state. These data ascertain two points: (1) there are three structural states with the cooling rates ranging from 5 1013 K/s to 5 1011 K/s; and (2) the smaller cooling rate leads to the more stable structure with larger degree of order. In order to investigate the size effect on the resultant solidification microstructures, we simulated the cooling process of a lager system with 48,668 atoms. The results of the larger system are almost the same as the system of 4000 atoms (not shown here), which indicates that the simulated system with 4000 atoms under periodic boundary conditions in this paper is large enough to investigate the thermodynamics and structural properties. During the process of glass formation, the atomic potential energy and volume decrease almost linearly with decreasing temperature except the temperature region near the glass transition temperature (Tg), whilst the slopes of the potential energy (or volume) versus temperature are different for the high temperature region above Tg and the low temperature region below Tg. According to this liquid-to-glass transformation character, the glass transition temperature Tg can be determined with the inflexion occurred in the potential energy (or volume) versus temperature curve. As such, the Tg corresponding to different cooling rates was estimated: Tg1 ¼ 570 K for Q1, Tg2 ¼ 525 K for Q2 and Tg3 ¼ 520 K for Q3, identifying the fact that the higher cooling rate leads to the higher Tg [18]. The higher cooling rate results in shorter time for atomic relaxation, thus leading to formation of the glass at a higher temperature than the lower cooling rate. During the crystallization process, there is a sudden drop of the potential energy and volume near the temperature of 950 K, which is caused by the first-order transition of crystallization. The starting temperature of potential energy decrease is the onset temperature of crystallization. From Fig. 1, it is noticed that the drop of potential energy and volume is drastic for Q9 and Q8 while gradual for Q7. In addition, below 500 K their energy and volume have almost the same values and decrease linearly with temperature for Q7, Q8, and Q9. Careful analyses reveal that the temperature regions for the crystallization process corresponding to different cooling rates are 950w900 K for Q9, 950w850 K for Q8, and 900w500 K for Q7, indicating that the temperature region increases with increasing cooling rate. The resultant structure during liquid-solid phase transformation is controlled by the cooling rate to a large extent. At slow cooling rates the atoms in liquids have enough time to relax to decrease free energy of the system as far as possible, and then the nucleation and growth occurs, resulting in the formation of crystalline solids. As the cooling rates increase, the relaxation time allowed for atoms becomes shorter and shorter, and thus the incompletely-relaxed atomic configuration is eventually frozen, leading to formation of the amorphous state. From above results, one can determine that the critical cooling rate for glass-forming is about 1.5 1013 K/s whilst the critical cooling rate for crystallization is about 4 1012 K/s. As the cooling rates range from 6 1012 K/s to 1 1013 K/s, one can obtain the transitional structure state between amorphous and crystal states. 3.2. Pair distribution function In the MD simulation, the PDF is a very important tool to probe the structural changes of a system from the statistical point of view. For a monatomic system, the PDF can be defined as [28e30],
Fig. 2. Pair distribution function (PDF) of molecular dynamics (MD) simulation during different cooling rates of Ni simulated system at 300 K. (a) All cooling rates. (b) Cooling rates at 5 1013 K/s, 1 1013 K/s and 5 1011 K/s. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).
V gðrÞ ¼ 2 N
n X nðrÞ p 4 r 2 Dr i¼1
! (2)
where V is the volume of simulation cell, N denotes the number of atoms, and n(r) is the number of atoms around a central atom within the distance interval between r and r þ Δr .The PDF denotes the probability of finding another atom at the radius of r around a given center atom.
Table 1 The position (r21, r22) and intensity (g21, g22) of the two subpeaks for the splitting second peak on pair distribution function under different cooling rates. Cooling rate (K/s) 13
5 10 2 1013 1.5 1013 1 1013 8 1012 6 1012 4 1012 1 1012 5 1011
A) r21 (
g21
r22 ( A)
g22
g21/g22
g21/g22
4.35 4.36 4.36 4.37 4.37 4.37 4.37 4.37 4.38
1.55 1.63 1.67 2.28 2.17 2.09 3.70 3.94 3.63
4.83 4.86 4.86 4.94 4.91 4.91 5.03 5.03 5.03
1.38 1.34 1.33 1.22 1.22 1.19 1.40 1.44 1.37
1.12 1.22 1.26 1.87 1.78 1.76 2.64 2.74 2.65
1.2 0.1
1.8 0.1
2.7 0.1
F. Li et al. / Intermetallics 19 (2011) 630e635
Fig. 2a displays the PDFs of Ni simulated system at 300 K obtained from different quenching rates. One can see that the first peak of PDF becomes higher and narrower while the first valley becomes deeper as the cooling rate decreases. The second peak of the PDF appears to split into two subpeaks for the higher cooling rates (Q1, Q2 and Q3), which indicates that the non-equilibrium phase obtained by fast quenching is amorphous state [31e33]. Also, the right subpeak intensity of the splitting peak (g22) is somewhat
633
lower than the left one (g21). As the cooling rate decreases to Q4 ¼ 1 1013 K/s, a new peak (shoulder peak) locating at w3.61 Å emerges as shown in Fig. 2b, which means that the ordered structure has been formed [34]. This shoulder peak becomes intensified and other new peaks subsequently emerge at larger r values in the PDF curves for the lower cooling rates (Q7, Q8 and Q9). These peaks correspond to the PDF of typical fcc crystal structure (Fig. 2a), indicating that the 3D periodic fcc structure has formed.
Fig. 3. Molecular dynamics (MD) simulated instantaneous atomic images of the Ni system at 300 K resulted from different cooling rates. (a) Amorphous structure at Q3 ¼ 1.5 1013 K/s. (b) 1D periodic structure at Q4 ¼ 1 1013 K/s, the green arrow denotes the direction of the 1D periodic structure. (c) Detailed atomic configuration of the selected plane (blue atoms) in (b). (d) 2D periodic structure at Q6 ¼ 6 1012 K/s, the green and blue arrows indicate the two directions of the 2D periodic structure. (e) 3D crystal structure at Q9 ¼ 5 1011 K/s. All the insets in the top right corner of these figures display the corresponding fast Fourier transformation (FFT) patterns of the atomic configurations for different cooling rates. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).
634
F. Li et al. / Intermetallics 19 (2011) 630e635
Specifically, for the lowest cooling rate Q9 (5 1011 K/s), the positions and even the relative values of all peaks are almost the same as those in the fcc crystal, implying the resultant solidification structure is close to a perfect fcc crystal. From Fig. 2a, one can see that the splitting second peak in PDF of amorphous state gradually decomposes into two separate peaks of crystal with decreasing cooling rate. The red solid and black dotted lines show the evolution trend. The positions (r21, r22) and intensities (g21, g22) of these two subpeaks which correspond separately to the third and the fourth peaks in the crystal structure are listed in Table 1. The results reveal that the first subpeak position (r21) and the second subpeak position (r22) shift toward the larger distances with decreasing cooling rate. Interestingly, the intensity ratios of these two peaks (g21/g22) for these three structure states correspond to three characteristic values 1.2, 1.8 and 2.7, respectively. Meanwhile, the larger characteristic value corresponds to the more ordering structure. The result suggests that the ratio of g21 to g22 could act as an indicator to determine the structure state of a material. 3.3. Configuration images analysis In order to obtain a global picture of the structural characteristics in the 3D space, the atomic configuration images resulted from different cooling rates were viewed along different directions by rotating. Fig. 3a shows the atomic configuration of Q3 at ambient temperature. It exhibits the typical amorphous structure characteristics without long-range order, which is consistent with the result of PDF. Moreover, we use fast Fourier transformation (FFT) patterns to characterize the structural features of the simulated system and the FFT formula is expressed as
ZN XðkÞ ¼
xðrÞexpði2pkrÞdr
(3)
N
where x(r) and X(k) separately denote the images in real space and reciprocal space. Parameters r and k represent vectors in real space and reciprocal space, respectively. The FFT pattern at the top right of Fig. 3a also verifies the long-range disordered feature of the atomic configuration for Q3. With the cooling rate decreasing to Q4, it is found that there are a series of parallel planes along one exclusive direction (its periodic direction as illustrated by the green arrow in Fig. 3b). This atomic configuration is similar with our previous finding in studying the crystallization of amorphous alloys [35,36], where this structure was referred to as one-dimensional (1D) periodic structure. The periodicity of this 1D structure can be characterized by two symmetrical diffraction spots (the inset at the top right of Fig. 3b) and the line connecting the two spots is perpendicular to the close-packed planes. One plane (the blue plane in Fig. 3b) was picked up from the 1D configuration to investigate the atomic packing on these parallel planes, and its detailed atomic distribution is shown in Fig. 3c. It can be seen that the atom arrangement exhibits a typical (111) close-packed plane in the fcc lattice, which is consistent with our previous work [35,36]. When the cooling rates further decrease to Q6, two-dimensional (2D) periodic structure is formed as demonstrated in Fig. 3d. In the 2D structure, the angle between the two periodic directions was measured to be about 60 (the periodic directions as displayed by the blue and green arrows in Fig. 3d). The FFT image of this ordered assemblage shows a typical diffraction pattern of 2D periodicity, as shown in the upper corner of Fig. 3d. The diagonal lines connecting these four diffraction spots are perpendicular to the two closepacked planes and the angle between the two diagonal lines was also about 60 . It should be noted that 1D and 2D periodic
structures are local and cannot run through the entire simulated Ni model. As for the atomic configurations resulted from the lowest cooling rate Q9, the 3D periodic crystal structure can be obviously observed (Fig. 3e). The FFT pattern also identifies the 3D periodicity (as shown at the top right of Fig. 3e). The lines connecting two symmetrical diffraction spots are parallel to the periodicity direction of the 3D periodicity structure, respectively. The above configurational analyses results are consistent with the analysis of potential energy, molar volume and PDF. 3.4. Bond pair analysis Honeycutt and Anderson (HA) bond-type index method [21,37,38] is capable to describe and discern the concrete relationship of an atom with its near neighbors. According to their formula, the local environment can be represented by the indices i, j, k and l. When the distance between atoms A and B is less than their “bonding” distance, i ¼ 1, otherwise i ¼ 2; j denotes the atom number of forming bonds with both atom A and atom B; k is the bond number of forming bonds among the above j atoms; and l is a parameter used to distinguish local structures when k is the same. Different structures are characterized by different representative HA bond-types: 1551, 1541, and 1431 pairs corresponding to icosahedral ordering are characteristic of typical liquid and amorphous states; the 1421 pair represents the fcc local structure; 1441 and 1661 pairs are characteristic of body-centered cubic (bcc) crystal; and the 1422 pair characterizes hexagonal close-packed (hcp) crystal. Fig. 4 represents the relative content of different kinds of bond pairs for the atomic configurations solidified at various cooling rates. The total concentration of these seven bond-types is over 88% for all the cooling rates, and there is an obvious dependence of the fraction of bond-types with the cooling rate. It is seen that the content of 1551, 1541 and 1431 bond pairs which are associated with the icosahedral ordering for runs Q1, Q2 and Q3, is remarkably greater than those of the others. The total amount of 1551, 1431 and 1541 type bond pairs sharply decrease from more than 60% to about 1% when the cooling rates fall from Q1 to Q9, implying that the icosahedral bond pairs are dominant in the amorphous structure. It is noted that the content of 1431 and 1541 bond pairs is larger than that of 1551 bond pair, indicating that there are more defected icosahedra than ideal icosahedra in the present amorphous Ni. However, the fraction of 1551 bond pair in transitional structures (Q4, Q5, and Q6) is as small as 2% while the total content of 1431 and 1541 pairs is more than 16%. Meanwhile, the content of 1421 bond pair is over 81% when the cooling rates drop to Q7, Q8 and Q9,
Fig. 4. Variation of the fractions of HoneycutteAnderson (HA) indices of Ni simulated system at 300 K during different cooling rates. The insets show the typical icosahedron and fcc structures. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).
F. Li et al. / Intermetallics 19 (2011) 630e635
indicating the structure basically becomes typical fcc crystal structure. As for the 1441 and 1661 bond pairs, their total content is less than 5% for all the cooling rates. It can be seen that the crystallization degree of a system can be described by the HA bond-type index method more clearly than that described by the PDF (Fig. 2a). In our study, the fraction of the 1422 bond-type does not monotonously increase with decreasing the cooling rate, which is coincident with the report in literature [19]. The concentration of the 1422 bond pair in the 1D and 2D transitional structures (Q4, Q5 and Q6) is about 20%, which is larger than that in the amorphous and crystal states. Moreover, the fraction of 1421 bond pair is about 45% in 1D and 2D periodic structures. In terms of the bond pair analyses and the instantaneous configuration image of a close-packed plane (Fig. 3c), we can reasonably infer that the fcc- and hcp-type closepacked planes make up the 1D and 2D transitional structures. 4. Conclusions The present MD simulations demonstrate that the cooling rate plays a critical role in final solidification microstructures of the liquid Ni. With the cooling rates decreasing from 5 1013 K/s to 5 1011 K/s, the final atomic structures of Ni at ambient temperature evolve from the amorphous to the fcc crystal states mediated by transitional structures with the characteristics of 1D and 2D periodicities. Further studies reveal that fcc- and hcp-type close-packed planes constitute the 1D and 2D periodic structures. The critical cooling rate to form the amorphous state is determined to be about 1.5 1013 K/s, whilst the critical cooling rate for crystallization is about 4 1012 K/s. In addition, the two subpeaks of the splitting second peak on the PDF gradually decompose into two separate peaks when the final solidification structure transforms from the amorphous to crystalline states with decreasing cooling rates. The intensity ratios of these two peaks g21/g22 for the amorphous, the transitional structure and the crystalline states correspond to three different characteristic values, respectively. Therefore, the value of g21/g22 could act as a parameter to characterize the structure state. Acknowledgements This work is supported by the grants from National Basic Research Program of China (No.: 2011CB605504), National Natural Science Foundation of China (Nos.: 50871054 and 50901006), Research Fund for the Doctoral Program of Higher Education of
635
China (Nos.: 20093219110035 and 20090006120025), Program of Introducing Talents of Discipline to Universities (Project No. B07003).
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38]
Mandell MJ, Mctague JP. J Chem Phys 1976;64:3699e702. Lu J, Qi DW. Phys Lett A 1991;157:283e5. Li M, Johnson WL. Phys Rev Lett 1993;70:1120e3. Coluzzi B, Parisi G, Verrocchio P. Phys Rev Lett 2000;84:306e9. Debenedetti PG, Stillinger FH. Nature 2001;410:259e67. Chan WL, Averback RS, Cahill DG, Ashkenazy Y. Phys Rev Lett 2009;102: 095701. Inoue A, Zhang T, Masumoto T. Mater Trans JIM 1990;31:425e8. Greer AL. Science 1995;267:1947e53. Johnson WL. Mater Res Bull 1999;24:42e56. Löffler JF. Intermetallics 2003;11:529e40. Wang WH, Dong C, Shek CH. Mater Sci Eng R Rep 2004;44:45e89. Hng HH, Li Y, Ng SC, Ong CK. J Non-Cryst Solids 1996;208:127e38. Takeuchi A, Inoue A. Mater Sci Eng A 2001;304e306:446e51. Hildal K, Sekido N, Perepezko JH. Intermetallics 2006;14:898e902. Li Q. Mater Lett 2007;61:3323e8. Huang YJ, Shen J, Chen JJJ, Sun JF. J Alloys Compd 2009;477:920e4. Hsu CS, Rahman A. J Chem Phys 1979;70:5234e40. ın T, Kimura Y, Goddard WA. Phys Rev B 1999;59:3527e33. Qi Y, Çag Liu CS, Xia JC, Zhu ZG, Sun DY. J Chem Phys 2001;114:7506e12. Sheng HW, He JH, Ma E. Phys Rev B 2002;65:184203. Sheng HW, Luo WK, Alamgir FM, Bai JM, Ma E. Nature 2006;439:419e25. Kazanc S. Comput Mater Sci 2006;38:405e9. Qi L, Dong LF, Zhang SL, Cui ZQ, Ma MZ, Jing Q, et al. Comput Mater Sci 2008;42:713e6. Fang HZ, Hui X, Chen GL, Öttking R, Liu YH, Schaefer JA, et al. Comput Mater Sci 2008;43:1123e9. Cleri F, Rosato V. Phys Rev B 1993;48:22e33. Rifkin J. XMD-molecular dynamics program, version 2.5.32; 2002. Willam WD. Materials science and engineering: an introduction. New York: John Wiley & Sons, Inc.; 1994. Allen MP, Tildesley DJ. Computer simulation of liquids. Oxford: Oxford University Press; 1987. Waseda Y. The structure of non-crystalline materials. New York: McGraw-Hill; 1980. Jin W, Kalia RK, Vashishta P, Rino JP. Phys Rev B 1994;50:118e31. Finney JL. Nature 1977;266:309e14. Sachdev S, Nelson DR. Phys Rev Lett 1984;53:1947e50. Zallen R. The physics of amorphous solids. New York: Wiley-Interscience; 1983. Li H. Phys Rev B 2003;68:024210. Liu XJ, Chen GL, Hui XD, Hou HY, Yao KF, Liu CT. J Appl Phys 2007;102:063515. Liu XJ, Chen GL, Hou HY, Hui XD, Yao KF, Lu ZP, et al. Acta Mater 2008;56: 2760e9. Honeycut JD, Anderson HC. J Phys Chem 1987;91:4950e63. Lo YC, Huang C, Ju SP, Du XH. Phys Rev B 2007;76:024103.