Temperature effects on atomic pair distribution functions of melts J. Ding, M. Xu, P. F. Guan, S. W. Deng, Y. Q. Cheng, and E. Ma Citation: The Journal of Chemical Physics 140, 064501 (2014); doi: 10.1063/1.4864106 View online: http://dx.doi.org/10.1063/1.4864106 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/140/6?ver=pdfcov Published by the AIP Publishing
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 66.173.107.86 On: Wed, 12 Mar 2014 20:35:03
THE JOURNAL OF CHEMICAL PHYSICS 140, 064501 (2014)
Temperature effects on atomic pair distribution functions of melts J. Ding,1,a) M. Xu,2 P. F. Guan,1,3 S. W. Deng,1,4 Y. Q. Cheng,5 and E. Ma1 1
Department of Materials Science and Engineering, Johns Hopkins University, Baltimore, Maryland 21218, USA 2 I. Physikalisches Institut IA, RWTH Aachen University, Aachen 52056, Germany 3 Beijing Computational Science Research Center, Beijing 100086, China 4 Department of Chemistry, East China University of Science and Technology, Shanghai 200237, China 5 Chemical and Engineering Materials Division, Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, USA
(Received 27 November 2013; accepted 17 January 2014; published online 10 February 2014) Using molecular dynamics simulations, we investigate the temperature-dependent evolution of the first peak position/shape in pair distribution functions of liquids. For metallic liquids, the peak skews towards the left (shorter distance side) with increasing temperature, similar to the previously reported anomalous peak shift. Making use of constant-volume simulations in the absence of thermal expansion and change in inherent structure, we demonstrate that the apparent shift of the peak maximum can be a result of the asymmetric shape of the peak, as the asymmetry increases with temperature-induced spreading of neighboring atoms to shorter and longer distances due to the anharmonic nature of the interatomic interaction potential. These findings shed light on the first-shell expansion/contraction paradox for metallic liquids, aside from possible changes in local topological or chemical short-range ordering. The melts of covalent materials are found to exhibit an opposite trend of peak shift, which is attributed to an effect of the directionality of the interatomic bonds. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4864106] I. INTRODUCTION
The liquid structure is key to understanding melting, glass transition, crystallization, and many other fundamental problems in materials science.1–4 Thanks to recent progress in experimental and computer simulation techniques, the atomic structure of many liquids have been characterized and correlated with macroscopic properties.5–11 In these studies, atomic pair distribution functions (PDFs) have been widely used to describe the amorphous structure.11 An intriguing behavior noticed in some of these studies is an anomalous temperature dependence of the first peak position in the PDF of liquids: the position (distance r) corresponding to the maximum intensity of the first peak shifts towards shorter r (hereafter abbreviated as “shift of the first peak”) with increasing temperature. This skewing of the first peak has been observed in a number of metallic melts, such as single-component Cu, Al, Zn, Ni, Au liquids12–16 and alloy liquids.17–19 Liquids simulated using the classical Lennard-Jones potential also follow this trend.20 In contrast, in liquids of covalent materials such as Si, Ge, Ge1 Sb2 Te4 , and polymers,21–26 the first peak shifts in the opposite direction with elevating temperature. While the latter behavior is expected from thermal expansion, that of the metallic liquids is counter-intuitive and therefore surprising. A recent study combining x-ray diffraction and molecular dynamics simulation attributed the anomalous temperature dependence of the first peak to a “negative expansion” of interatomic distances in metallic melts.12 Specifically, they argued that low-coordinated atomic clusters become more preferable at higher temperatures, and this change in local topological a) [email protected]
0021-9606/2014/140(6)/064501/8/$30.00
order leads to a contraction of the average atomic distance between the central atom and the surrounding nearest neighbors. In this paper we analyze several typical liquids, to identify potential factors responsible for the shift of the first PDF peak with increasing temperature. These liquids are studied using classical molecular dynamics (MD) simulations27, 28 or ab initio molecular dynamics29 calculations. Detailed structural analysis will be conducted, in light of the potential energy landscape (PEL) picture,30, 31 to gain insight into the temperature dependence of the atomic configurations and pair distributions in liquids. The questions to be addressed include: (i) Is the r position of the first peak maximum in PDFs of liquids a suitable measure to gauge the changes in interatomic distances and in the inherent atomic configuration? (ii) Would the coordination numbers (CNs) actually decrease with increasing temperature in metallic liquids? (iii) Can the peak shift “anomaly” in metallic liquids originate from the shape evolution of the peak, and particularly its increasing degree of asymmetry associated with anharmonicity at elevated temperatures? (iv) And why liquids of covalent materials appear “normal,” in terms of their T-dependence of the first peak position? II. METHODS
Molecular dynamics simulations of several representative liquids were carried out. The metallic liquids studied include Cu, Al, Cu64 Zr36 , employing embedded atom method (EAM) potentials.32, 33 Covalent Si and Ge2 Sb2 Te5 (GST) were also melted at high temperatures, to compare the other liquids. The
140, 064501-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 66.173.107.86 On: Wed, 12 Mar 2014 20:35:03
064501-2
Ding et al.
J. Chem. Phys. 140, 064501 (2014)
Stillinger-Weber potential22 was used for Si, which includes both a pair interaction term, U2 (r), and a three-body interaction term, U3 (ri , rj , rk ). The reduced pair interaction is A(Br −p − r −q ) exp[(r − a)−1 ] r < a, U2 (r) = (1) 0 r≥a and the three-body term (i.e., the angle-dependent term) can be expressed as U3 (ri , rj , rk ) = h(rij , rik, θj ik ) +h(rj i , rj k, θij k ) + h(rki , rkj, θikj ),
(2)
where h(rij , rik, θj ik ) = λ exp[γ (rij − a)−1 + γ (rik − a)−1 ] × (cos θj ik + 13 )2 . Details of each parameter can be found in Ref. 12. We also utilized the pair interaction in the StillingerWeber potential (i.e., without the three-body term) to simulate a virtual “pair-potential liquid Si” (ppl-Si) to assess the influence of the bond directionality. Periodic boundary conditions (PBCs) were applied in all three directions. Cu, Al, Cu-Zr, Si, and ppl-Si liquids each contains 128 000 atoms; the system sizes are sufficiently large to minimize the artifacts associated with the PBC. Two different ensembles (NPT and NVT) were applied to prepare liquid configurations at each temperature of interest. The liquid samples were equilibrated at near or above the melting temperature, under Nose-Hoover thermostat (for NPT and NVT) and barostat (for NPT) with zero external pressure.28 The density for each sample under NVT ensemble was kept constant (the zero pressure density at the melting temperature). For each ensemble, 50 independent instantaneous liquid configurations (ILCs) were recorded and averaged. The corresponding inherent structures (ISs) were obtained by energy minimization using the conjugate gradient method. The GST liquid configurations were obtained from ab initio molecular dynamics simulations using the Vienna ab initio simulation package (VASP),29 with canonical (NVT) ensembles and temperature controlled by the Nose thermostat. Projector augmented-wave (PAW) method with the Perdew-Wang exchange-correlation potentials were adopted.34, 35 The simulation was performed at the point only. The GST samples each had 378 atoms (84 Ge atoms, 84 Sb atoms, and 210 Te atoms under PBCs) and were melted at 3000 K for 2000 time steps (6 ps) and then quenched to a specific temperature of interest, followed by relaxation. Six independent configurations were then selected for structural analysis. To investigate the distance-dependent dilation strain between the instantaneous liquid configurations and the corresponding inherent structure, we define σ thermal (r) following.36–38 By grouping the atomic pairs within the shell between (r − r) and (r + r) in the instantaneous liquid configurations, we solve the affine transformation matrix Jr which best map the following:
Jr rjIiLC → rjIiS ; r − r < rj0i ≤ r + r.
FIG. 1. PDFs of Al liquids at several temperatures under zero external pressure (NPT ensemble). The inset shows a magnified view of their first peaks.
yy
ηxx +η +ηzz
σthermal (r) = r 3r r . In our definition, positive σ thermal (r) means that the shell at the distance of r in the ILCs will dilate upon energy minimization to reach the corresponding IS. III. RESULTS AND DISCUSSIONS
First, we examine the temperature dependence of PDFs in metallic melts. In Fig. 1(a), PDFs of ILCs for Al melts from 900 K to 1900 K under zero external pressure (NPT ensemble) are plotted. The first peak, blown-up in the inset, exhibits the “anomalous” behavior noticed in a number of previous experiments and simulations:12–16 the position corresponding to the maximum peak height shifts to shorter distance with increasing temperature. We also simulated liquids with a fixed density. Figure 2 presents the PDFs of both ILC and IS of Al melts under the constant-volume (NVT ensemble) simulation. As expected, the PDFs of the ISs at various temperatures all collapse onto one single profile, exhibiting
(3)
The Lagrangian strain matrix can then be calculated as ηr =
1 Jr JrT − I . 2
(4)
FIG. 2. PDFs of Al liquids and their inherent structures at various temperatures (simulated under NVT ensemble, i.e., the density is fixed and temperature independent).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 66.173.107.86 On: Wed, 12 Mar 2014 20:35:03
064501-3
Ding et al.
FIG. 3. Schematic showing instantaneous liquid configurations (solid circles) and the corresponding inherent structures (open circles), from the perspective of potential energy landscape (PEL).
negligible temperature dependence. This is consistent with the PEL theory, i.e., equilibrium liquids at a constant density are iso-configurational.22, 30, 31 With no volume expansion nor evolution of inherent structure, here the first peak position in PDFs of ILCs clearly shifts to the left with increasing temperature (similar to the NPT results in Fig. 1 and Ref. 22). As such, Fig. 2 demonstrates that the temperature-dependent shifts in the PDFs of the ILC can arise from thermal effects alone, in the absence of any appreciable evolution of the IS22, 30 or positive/negative thermal expansion.
J. Chem. Phys. 140, 064501 (2014)
To illustrate the difference between the ILC and its corresponding IS state, Fig. 3 schematically depicts a PEL,30, 31 with solid and open circles representing the ILCs and ISs, respectively. The IS is the local minimum of the potential energy basin in which the ILC resides, after removing the temperature effects through energy minimization (the IS may be regarded as the go-to configuration at zero temperature). From this perspective, the results in Fig. 1 indicate that it is the increased thermal excitation at higher temperatures that induces a re-distribution of atoms over a wider range of distances to cause an anomalous shift of the first peak maximum towards the shorter r side. This latter conjecture can be shown to hold, in the absence of changes in system IS/density, as done in Fig. 2. Not all liquids behave this way, though. In Fig. 4, the PDFs of ILC and the corresponding IS for (a) Cu, (b) Cu64 Zr36 , (c) Si, and (d) Ge2 Sb2 Te5 39 are presented. For Cu and Cu64 Zr36 ILCs, their first peak positions in PDFs shift to shorter r with increasing temperature, just like Al (see Figs. 1 and 2) and other metallic liquids. On the contrary, the first peak position in PDFs of Si and Ge2 Sb2 Te5 liquids moves to larger r with rising temperature, which agrees with other experimental and simulation data of molten covalent Si, Ge, and polymers.22–26 Meanwhile, the PDFs of IS for each of the liquid systems in Fig. 4 always overlap (there are
FIG. 4. PDFs of (a) Cu, (b) Cu64 Zr36 , (c) Si, and (d) GST liquids and their corresponding inherent structures at various temperatures at constant density (NVT ensemble).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 66.173.107.86 On: Wed, 12 Mar 2014 20:35:03
064501-4
Ding et al.
minor variations for GST liquids due to the limited box-size and number of configurations in DFT calculation). Also, by comparing the PDFs of ILC and IS in Fig. 4 (and Fig. 2), it can be seen that for metallic liquids, the first peak position in PDFs of ILC is on the left-hand (lower r) side relative to the peak position in the IS PDFs (note that the splitting of the first peak in Cu-Zr inherent structure comes from the different atomic radius of Cu and Zr atoms), but for Si and GST liquids the peak position of the ILC PDF is on the right-hand (larger r) side. Since IS can be regarded as a configuration at zero temperature, the above is also consistent with the observation that with increasing temperature, the ILC peak position shifts towards the left for metallic liquids, but to the right for covalent liquids. In other words, the position corresponding to the IS corresponds to the trend extrapolated to zero temperature. In liquids, the nearest neighbors (NNs) and CNs are not unambiguously defined. In order to evaluate NN or CN, a distance cutoff at the minimum after the first peak in PDF is widely employed. But caution should be exercised when it is used to analyze subtle changes in NN or CN. The distance cutoff is not set at zero PDF intensity, and the minimum (the dip) in fact corresponds to a fairly high intensity (see Fig. 1), especially when considering that for high-temperature liquids the peak intensity distribution spans a rather wide r range. In our constant-density simulations, since ILC at all temperatures is iso-configurational with virtually identical IS, the CNs of their IS should also be the same no matter which method is employed to determine the NNs. In contrast, if we analyze ILC using the distance cutoff at first minimum, the CN will decrease with increasing temperature (not shown here), as observed in Ref. 12. This is because from IS to ILC, and from ILC at a lower temperature to that at a higher temperature, the profile (shape) of the PDF peak changes (broadens) due to thermal activation, such that the pre-set cutoff misses part of the first neighbors that have moved to longer r. Therefore, the apparent change in CN for the ILCs, when determined by the distance cutoff method, is merely a dynamical effect due to thermally induced spreading. The same is true for the “first peak position,” as will be explained next. In Fig. 4, it is observed that the profile of the first peak in PDFs becomes increasingly asymmetric with increasing temperature. To evaluate the asymmetry, we employ the ratio of rB /rA , with the two distances rB and rA defined at the full width of half maximum height of the first PDF peak, see inset in Fig. 5 for a Cu liquid at 1300 K. For a symmetric peak, rB /rA will be equal to 1. Fig. 5 shows that the rB /rA for Al, Cu, and Cu-Zr liquids (Si and GST liquids are excluded for this analysis because their first peaks are not sufficiently sharp) is always greater than 1 and increases with increasing temperature. This clearly indicates that higher temperature leads to higher asymmetry of the first peak in PDFs for equilibrium metallic liquids. Specifically, for the first PDF peak the curve at shorter distances (the left wing of the peak) is steeper than that at larger distances (right wing). This is due to the shape of the typical interatomic interaction potential for pair interactions (the potential energy versus r curve):20, 22 the anharmonicity of the potential grows as atoms reach higher energy
J. Chem. Phys. 140, 064501 (2014)
FIG. 5. Temperature dependence of the asymmetry of the first peak in PDFs (evaluated by rB /rA ) for Cu, Al, and Cu-Zr liquids (constant density). The inset schematically shows the definition of rB and rA (at the full width of half maximum of the first peak for Cu liquid at 1300 K; purple circles).
levels at elevated temperatures. More discussions on this point will be presented later. To see the re-distribution due to temperature, we compare the ILC and the IS, at a constant density, to observe how much the atoms move upon energy minimization. The σ thermal (r) defined in Sec. II can measure the distance-dependent dilation strain from ILC to the corresponding IS. For instance, a positive σ thermal (r) means that the shell at the distance of r in ILC will dilate when reaching the corresponding IS and vice versa. As presented in Fig. 6(a), σ thermal (r) of Cu liquids at 1300 K and 2500 K are very similar and the dilation strain ranges between −6% and 21%. In other words, some interatomic distances would dilate as much as 21% upon energy minimization, while some others shrink by a few percent. Also, a zero value of σ thermal (r) in Fig. 6(a) means that a balance is reached between dilation and contraction for that shell at the distance of r. For comparison, Fig. 6(b) shows the average r number of atomic pairs within the distance of r, N (r) = 0 4π r 2 ρ0 g(r)dr(where ρ 0 is the average number density) for Cu liquids at 1300 K, 2500 K, and IS (see left axis), as well as the PDF of IS for Cu liquid (see right axis). The dashed lines in Fig. 6 indicate the locations where σ thermal (r) is equal to zero. We observe that (i) σ thermal (r) = 0 corresponds to the crossing between the N(r) of IS and the N(r) at high temperatures and (ii) the zero value of σ thermal (r) also corresponds to either the peak or valley position (except R2 ) in the PDF of (r) < 0 corresponds IS. Specifically, σ thermal (r) = 0 and dσthermal dr (r) >0 to peak position, whereas σ thermal (r) = 0 and dσthermal dr coincides with the minimum valley in PDF in Fig. 6(b). We observe that the position of R2 does not directly point to peak position; this is due to the splitting of this peak. But R2 /R1 ≈ 1.8 is consistent with the prediction by spherical-periodic ,40 where n is an integer, order (SPO), RnSP O /R1SP O = 15 + 4n 5 as well as the average position of the second atomic shell as calculated for five different amorphous alloys.41 Therefore, R2 , as an average position of the second shell, is in line with the interpretation of R1 and R3 (R3 /R1 ≈ 2.6, also in agreement with SPO prediction). In sum, the observations above
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 66.173.107.86 On: Wed, 12 Mar 2014 20:35:03
064501-5
Ding et al.
FIG. 6. (a) Distance dependence of the dilation strain (σ thermal ) from instantaneous liquid configuration to the corresponding inherent structure in Cu liquids at 1300 K and 2500 K. The dashed lines indicate the position r with zero value of σ thermal . (b) N (r) = 0 4π r 2 ρ0 g(r)dr for Cu at 1300 K, 2500 K and the corresponding inherent structure (left axis), and the PDF of inherent structure of Cu liquids at 1300 K (right axis).
provide a general picture of the mapping between an ILC and its corresponding IS: (i) in ILC, the interatomic distances in IS spread out to become both shorter and longer under thermal agitation; (ii) the peaks (centers of atomic shells) and minima (valleys separating two neighboring shells) reach the specific condition of σ thermal (r) = 0 under thermal agitation. From the information above, two reasons come to mind that can explain why the position of the maximum intensity of the first peak in PDFs of metallic melts shifts to shorter r with increasing temperature. First, as shown in the preceding paragraph and Fig. 6, thermal agitation of atoms in metallic liquids leads to both shorter and longer bonds (i.e., negative and positive departure from the bond length in IS, RIS , which corresponds to the IS first peak position). The average bond length (Rave ), estimated considering CN to be the same as that in the corresponding IS, is plotted in Fig. 7(a) as a function of temperature for each of the two groups of atoms, one with shorter bonds (r < RIS ) and the other with longer ones (r > RIS ). Starting from the two data points at 0 K for IS, we observe that with increasing temperature, the two groups diverge, indicating an expanding atomic distance range. The Rave of the r > RIS group is larger than that of the r < RIS group. This, again, is because typical pair interactions are asymmetric and anharmonic, which becomes especially pronounced at high temperatures. The potential-distance curve is usually steeper on the shorter distance side but becomes much less so for longer r. In other words, the resistance the thermally excited atoms have to combat with is stiffer on the r < RIS side but lower on the r > RIS side. As a result, under thermal activation it is easier for atoms to move further away from the
J. Chem. Phys. 140, 064501 (2014)
central atom to reach r > RIS , with a distance range wider than that on the r < RIS side. To counter balance these latter atoms such that the system average is still at RIS , a higher density of atoms is expected to populate on the narrower r < RIS side. This asymmetric atomic distribution is illustrated in Fig. 7(b), reflected by a higher PDF peak intensity at smaller r in the left wing, together with broadening towards larger r on the right side. One can see this by hypothetically assuming the same number of atoms drifting to the left of RIS and to the right of RIS . Those on the right tend to take r positions further away from RIS , while those on the left reside closer to RIS over a smaller r range. Consequently the number density appears to be higher on the smaller r side. This is the reason for the enhanced asymmetry of the first PDF peak at high temperatures (also seen in Fig. 5). Such a skewing of the peak results in an apparent shift of the peak maximum position to lower r with increasing temperature, even when the system has the same average bond length (RIS ) in our constant-volume (and constant IS and CN) simulation. A similar trend is expected under constant pressure conditions. Second, another contributor to the change of the first peak position/shape is the normalization of 4π r2 in the calculation of g(r). As g(r) = 4πrS(r) 2 ρ dr , where S(r) is the average 0 number of atoms within the shell with atomic distance between (r − 12 dr) and (r + 12 dr), the g(r) peak position will be different from that in the radial distribution function, RDF (= 4π r2 g(r)), as seen from the difference between RPDF and RRDF in Fig. 7(c). This difference between g(r) and 4π r2 g(r) increases with rising temperature (see the plot of RRDF -RPDF in Fig. 7(d)). For a given change in S(r) due to thermal motion, if it happens on the longer r side it would drop the PDF intensity more, when compared with the effect of the same change in S(r) on shorter r side. This is because RDF at higher r is divided by a larger normalization factor 4π r2 . This causes a relative enhancement of the intensity of the peak on the leftwing side (smaller r), contributing to the peak asymmetry and the apparent position shift of the maximum intensity of the peak. The last question to address is the origin of the different behavior observed in melts of covalent Si and GST. In Fig. 8, we plot the distance-dependent N(r) for Cu and Si at different temperatures. The average first peak positions (RIS ) in PDFs are chosen as the reference point. A clear difference between Cu and Si can be observed at the position of RIS and in the region immediately below RIS : for Cu we see more atoms than IS, whereas for Si there are fewer atoms than IS. In other words, there appears to be an additional driving force in liquid Si to move a higher fraction of surrounding nearest neighbor atoms to distances a bit longer than RIS than in liquid Cu. This is the reason why the position of the maximum of the first peak in the PDFs of Si (and GST and other covalent liquids) shifts to longer r with increasing temperature, which contrasts with the case in metallic melts such as Cu, Al, and Cu-Zr. This difference is rooted in an extra contribution from the directional covalent bonding. To see this we compare the standard Stillinger-Weber potential22 and a modified version without the three-body (directional) term, as described in Sec. II. The potential with pair interaction only is plotted in
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 66.173.107.86 On: Wed, 12 Mar 2014 20:35:03
064501-6
Ding et al.
J. Chem. Phys. 140, 064501 (2014)
FIG. 7. (a) Temperature dependence of the average interatomic distance for atom pairs with r ≤ RIS (shorter distance) and r ≥ RIS (longer distance) in Cu melts; (b) schematic showing increasing asymmetry of the first PDF peak with increasing temperature; (c) pair distribution function (PDF) compared with radial distribution function (RDF) for Cu liquid at 1300 K, as well as the corresponding ISs; and (d) temperature dependence of the difference between the first peak positions in RDF and PDF (denoted as RRDF − RPDF ).
the inset of Fig. 9 (solid line). This pair-potential simulated liquid without the three-body angular term is called ppl-Si hereafter for convenience, and the temperature dependence of the PDF of ppl-Si in reduced units is shown in Fig. 9. In contrast with that in Si, the first peak now also shifts to shorter r with increasing temperature, resembling the behavior in metallic liquids. This result clearly demonstrates the
role of the directionality of the covalent bonds. Specifically, liquid Si originates from covalent Si. In the melt the bonding becomes increasingly more metallic, while retaining certain degree of covalency. The liquid structure exhibits a more pronounced signature of ∼60◦ bond angle, with increasing deviation from its equilibrium bond angle (corresponding to cos (θ ) = 1/3).42 As seen in Eq. (2), the interaction energy would be
FIG. 8. Distance-dependent N(r) for ILC and IS of (a) Cu and (b) Si liquids. Colored regions are for r ≤ RIS . Note that when the number of neighbors is counted up to RIS (corresponding to N(r) at r = RIS ), the number of neighbors in the ILC is larger than that in the IS for Cu, but the opposite is true for Si.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 66.173.107.86 On: Wed, 12 Mar 2014 20:35:03
064501-7
Ding et al.
FIG. 9. PDFs of ppl-Si (potential modified to eliminate the three-body angular term) at various temperatures (in reduced units). The inset compares the modified potential (solid line) with the original Stillinger-Weber potential of Si.
elevated by the angular term even for unchanged bond length. For a simplistic configuration of only three atoms (see Fig. 9), which form an equilateral triangle with bond length r, their reduced potential-distance curve is plotted as dotted line in the inset of Fig. 9 if the angular term in Eq. (2) is included. A comparison with the curve from pair potential only (solid line) indicates that the angular term in Stillinger-Weber potential results in a shift of the potential-distance curve to higher r and makes it shallower on the high r side when compared with the pair interaction case. This moves some nearest neighbors to slightly longer distances. In other words, thermal effect in Si liquids enhances the deviation of the bond angle, such that atomic bonds tend to stretch to longer lengths. As a result of these atoms accumulated at such distances, the peak maximum in PDFs shifts to the right (longer r) with increasing temperature, in contrasts with the behavior in metallic liquids.
IV. CONCLUSION
In summary, we have employed classical and ab initio molecular dynamics simulations to systematically study the effects of temperature/thermal excitation on the atomic pair distributions in several representative liquids (both metallic and covalent). A comparative analysis of the ILC and IS has been conducted. Our results show that: (1) The seemingly anomalous peak shift with increasing temperature observed for metallic liquids is not truly anomalous, but intrinsic, reproducible, and predictable, see (2) below. It can happen even when the liquid experiences no change in inherent topological12 or chemical (such as in alloys with two or more constituent elements) short-range ordering. (2) One should study the entire peak profile, rather than just the position corresponding to the maximum intensity. The latter itself is not a meaningful and reliable metric to monitor and understand the changes in CN, bond length (interatomic distance), and local order. The
J. Chem. Phys. 140, 064501 (2014)
PDF as measured includes both configurational (inherent structure, local minimum in PEL) and thermal (deviation from the PEL minimum due to thermal activation) contributions. One intrinsic reason for the position/shape evolution is the thermally induced asymmetric redistribution of neighboring atoms to shorter and longer distances due to the anharmonic nature of the interatomic interaction potential. With increasing temperature, the PDF peak (the distribution of the atoms in the first shell) becomes increasingly more asymmetric. It would be inadvisable to draw conclusions based solely on a comparison of the positions of the intensity maximum of two peaks with different degree of asymmetry. (3) The opposite trend of first peak shift for metallic and covalent systems is systematic (although exceptions may exist, such as multi-component liquids or supercooled liquids under structural ordering). In particular, we note that the directionality term in the interaction potential influences the population and spreading of atomic distributions under thermal agitation, hence the direction and degree of peak shift with increasing temperature.
ACKNOWLEDGMENTS
J.D. and E.M. were supported at JHU by U.S.-DOEBES, Division of Materials Sciences and Engineering, under Contract No. DE-FG02-09ER46056. The computer simulations were performed using the NERSC supercomputers. Y.Q.C. was supported by the Scientific User Facilities Division, Office of Basic Energy Sciences, U.S. Department of Energy. 1 P.
G. Debenedetti and F. H. Stillinger, Nature (London) 410, 259 (2001). C. Frank, Proc. R. Soc. London, Ser. A 215, 43 (1952). 3 C. A. Angell, MRS Bull. 33, 544–555 (2008). 4 Z. H. Jin, P. Gumbsch, K. Lu, and E. Ma, Phys. Rev. Lett. 87, 055703 (2001). 5 Y. T. Shen, T. H. Kim, A. K. Gangopadhy, and K. F. Kelton, Phys. Rev. Lett. 102, 057801 (2009). 6 D. Holland-Moritz, D. M. Herlach, and K. Urban, Phys. Rev. Lett. 71, 1196 (1993). 7 J. C. Bendert, A. K. Gangopadhyay, N. A. Mauro, and K. F. Kelton, Phys. Rev. Lett. 109, 185901 (2012). 8 C. Tang and P. Harrowell, Nat. Mater. 12, 507–511 (2013). 9 J. Ding, Y. Q. Cheng, H. W. Sheng, and E. Ma, Phys. Rev. B 85, 060201 (2012). 10 J. Ding, Y. Q. Cheng, and E. Ma, Acta Mater. 61, 3130–3140 (2013). 11 Y. Q. Cheng and E. Ma, Prog. Mater. Sci. 56, 379–473 (2011). 12 H. B. Lou, X. D. Wang, Q. P. Cao, D. X. Zhang, J. Zhang, T. D. Hu, H. K. Mao, and J. Z. Jiang, Proc. Natl. Acad. Sci. U.S.A. 110, 10068–10072 (2013). 13 G. W. Lee, A. K. Gangopadhyay, K. F. Kelton, R. W. Hyers, T. J. Rathz, J. R. Rogers, and D. S. Robinson, Phys. Rev. Lett. 93, 037802 (2004). 14 P. Ganesh and M. Widom, Phys. Rev. B 74, 134205 (2006). 15 T. H. Kim and K. F. Kelton, J. Chem. Phys. 126, 054513 (2007). 16 N. A. Mauro, J. C. Bendert, A. J. Vogt, J. M. Gewin, and K. F. Kelton, J. Chem. Phys. 135, 044502 (2011). 17 N. A. Mauro, V. Wessels, J. C. Bendert, S. Klein, A. K. Gangopadhyay, M. J. Kramer, S. G. Hao, G. E. Rustan, A. Kreyssig, A. I. Goldman, and K. F. Kelton, Phys. Rev. B 83, 184109 (2011). 18 N. A. Mauro, W. Fu, J. C. Bendert, Y. Q. Cheng, E. Ma, and K. F. Kelton, J. Chem. Phys. 137, 044501 (2012). 19 K. Georgarakis, D. V. Louzguine-Luzgin, J. Antonowicz, G. Vaughan, A. R. Yavari, T. Egami, and A. Inoue, Acta Mater. 59, 708–716 (2011). 20 W. Kob and H. Andersen, Phys. Rev. E 51, 4626 (1995). 2 F.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 66.173.107.86 On: Wed, 12 Mar 2014 20:35:03
064501-8 21 T.
Ding et al.
H. Kim, G. W. Lee, B. Sieve, A. K. Gangopadhyay, R. W. Hyers, T. J. Rathz, J. R. Rogers, D. S. Robinson, K. F. Kelton, and A. I. Goldman, Phys. Rev. Lett. 95, 085501 (2005). 22 F. H. Stillinger and T. A. Weber, Phys. Rev. B 31, 5262 (1985). 23 S. Sastry and C. A. Angell, Nat. Mater. 2, 739–743 (2003). 24 A. Filipponi and A. D. Cicco, Phys. Rev. B 51, 12322 (1995). 25 K. S. Chang, Y. L. Wang, C. H. Kang, H. J. Wei, Y. H. Weng, and K. L. Tung, J. Membr. Sci. 382, 30 (2011). 26 O. Alexiadis and V. G. Mavrantzas, Macromolecules 46, 2450–2467 (2013). 27 D. C. Rapaport, The Art of Molecular Dynamics Simulation (Cambridge University Press, 2004). 28 M. P. Allen and D. J. Tidesley, Computer Simulation of Liquids (Clarendon Press, Oxford, 1989). 29 G. Kresse and J. Furthmüller, J. Comput. Mater. Sci. 6, 15–50 (1996). 30 F. H. Stillinger, Science 267, 1935–1939 (1995).
J. Chem. Phys. 140, 064501 (2014) 31 D.
Wales, Energy Landscapes: Applications to Clusters, Biomolecules and Glasses, 1st ed. (Cambridge University Press, 2004). 32 H. W. Sheng, M. J. Kramer, A. Cadien, T. Fujita, and M. W. Chen, Phys. Rev. B 83, 134118 (2011). 33 Y. Q. Cheng, E. Ma, and H. W. Sheng, Phys. Rev. Lett. 102, 245501 (2009). 34 P. E. Blöchl, Phys. Rev. B 50, 17953–17979 (1994). 35 G. Kresse and D. Joubert, Phys. Rev. B 59, 1758–1775 (1999). 36 M. L. Falk and J. S. Langer, Phys. Rev. E 57, 7192 (1998). 37 J. Li, Modell. Simul. Mater. Sci. Eng. 11, 173 (2003). 38 J. Ding, Y. Q. Cheng, and E. Ma, Appl. Phys. Lett. 101, 121917 (2012). 39 M. Xu, Y. Q. Cheng, H. W. Sheng, and E. Ma, Phys. Rev. Lett. 103, 195502 (2009). 40 P. Häussler, Phys. Rep. 222, 65 (1992). 41 J. Ding, Y. Q. Cheng, and E. Ma, “Atomic order beyond short-range scale in metallic glasses and liquids” (unpublished). 42 C. S. Liu, Z. G. Zhu, J. C. Xia, and D. Y. Sun, Phys. Rev. B 60, 3194 (1999).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 66.173.107.86 On: Wed, 12 Mar 2014 20:35:03
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 66.173.107.86 On: Wed, 12 Mar 2014 20:35:03
THE JOURNAL OF CHEMICAL PHYSICS 140, 064501 (2014)
Temperature effects on atomic pair distribution functions of melts J. Ding,1,a) M. Xu,2 P. F. Guan,1,3 S. W. Deng,1,4 Y. Q. Cheng,5 and E. Ma1 1
Department of Materials Science and Engineering, Johns Hopkins University, Baltimore, Maryland 21218, USA 2 I. Physikalisches Institut IA, RWTH Aachen University, Aachen 52056, Germany 3 Beijing Computational Science Research Center, Beijing 100086, China 4 Department of Chemistry, East China University of Science and Technology, Shanghai 200237, China 5 Chemical and Engineering Materials Division, Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831, USA
(Received 27 November 2013; accepted 17 January 2014; published online 10 February 2014) Using molecular dynamics simulations, we investigate the temperature-dependent evolution of the first peak position/shape in pair distribution functions of liquids. For metallic liquids, the peak skews towards the left (shorter distance side) with increasing temperature, similar to the previously reported anomalous peak shift. Making use of constant-volume simulations in the absence of thermal expansion and change in inherent structure, we demonstrate that the apparent shift of the peak maximum can be a result of the asymmetric shape of the peak, as the asymmetry increases with temperature-induced spreading of neighboring atoms to shorter and longer distances due to the anharmonic nature of the interatomic interaction potential. These findings shed light on the first-shell expansion/contraction paradox for metallic liquids, aside from possible changes in local topological or chemical short-range ordering. The melts of covalent materials are found to exhibit an opposite trend of peak shift, which is attributed to an effect of the directionality of the interatomic bonds. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4864106] I. INTRODUCTION
The liquid structure is key to understanding melting, glass transition, crystallization, and many other fundamental problems in materials science.1–4 Thanks to recent progress in experimental and computer simulation techniques, the atomic structure of many liquids have been characterized and correlated with macroscopic properties.5–11 In these studies, atomic pair distribution functions (PDFs) have been widely used to describe the amorphous structure.11 An intriguing behavior noticed in some of these studies is an anomalous temperature dependence of the first peak position in the PDF of liquids: the position (distance r) corresponding to the maximum intensity of the first peak shifts towards shorter r (hereafter abbreviated as “shift of the first peak”) with increasing temperature. This skewing of the first peak has been observed in a number of metallic melts, such as single-component Cu, Al, Zn, Ni, Au liquids12–16 and alloy liquids.17–19 Liquids simulated using the classical Lennard-Jones potential also follow this trend.20 In contrast, in liquids of covalent materials such as Si, Ge, Ge1 Sb2 Te4 , and polymers,21–26 the first peak shifts in the opposite direction with elevating temperature. While the latter behavior is expected from thermal expansion, that of the metallic liquids is counter-intuitive and therefore surprising. A recent study combining x-ray diffraction and molecular dynamics simulation attributed the anomalous temperature dependence of the first peak to a “negative expansion” of interatomic distances in metallic melts.12 Specifically, they argued that low-coordinated atomic clusters become more preferable at higher temperatures, and this change in local topological a) [email protected]
0021-9606/2014/140(6)/064501/8/$30.00
order leads to a contraction of the average atomic distance between the central atom and the surrounding nearest neighbors. In this paper we analyze several typical liquids, to identify potential factors responsible for the shift of the first PDF peak with increasing temperature. These liquids are studied using classical molecular dynamics (MD) simulations27, 28 or ab initio molecular dynamics29 calculations. Detailed structural analysis will be conducted, in light of the potential energy landscape (PEL) picture,30, 31 to gain insight into the temperature dependence of the atomic configurations and pair distributions in liquids. The questions to be addressed include: (i) Is the r position of the first peak maximum in PDFs of liquids a suitable measure to gauge the changes in interatomic distances and in the inherent atomic configuration? (ii) Would the coordination numbers (CNs) actually decrease with increasing temperature in metallic liquids? (iii) Can the peak shift “anomaly” in metallic liquids originate from the shape evolution of the peak, and particularly its increasing degree of asymmetry associated with anharmonicity at elevated temperatures? (iv) And why liquids of covalent materials appear “normal,” in terms of their T-dependence of the first peak position? II. METHODS
Molecular dynamics simulations of several representative liquids were carried out. The metallic liquids studied include Cu, Al, Cu64 Zr36 , employing embedded atom method (EAM) potentials.32, 33 Covalent Si and Ge2 Sb2 Te5 (GST) were also melted at high temperatures, to compare the other liquids. The
140, 064501-1
© 2014 AIP Publishing LLC
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 66.173.107.86 On: Wed, 12 Mar 2014 20:35:03
064501-2
Ding et al.
J. Chem. Phys. 140, 064501 (2014)
Stillinger-Weber potential22 was used for Si, which includes both a pair interaction term, U2 (r), and a three-body interaction term, U3 (ri , rj , rk ). The reduced pair interaction is A(Br −p − r −q ) exp[(r − a)−1 ] r < a, U2 (r) = (1) 0 r≥a and the three-body term (i.e., the angle-dependent term) can be expressed as U3 (ri , rj , rk ) = h(rij , rik, θj ik ) +h(rj i , rj k, θij k ) + h(rki , rkj, θikj ),
(2)
where h(rij , rik, θj ik ) = λ exp[γ (rij − a)−1 + γ (rik − a)−1 ] × (cos θj ik + 13 )2 . Details of each parameter can be found in Ref. 12. We also utilized the pair interaction in the StillingerWeber potential (i.e., without the three-body term) to simulate a virtual “pair-potential liquid Si” (ppl-Si) to assess the influence of the bond directionality. Periodic boundary conditions (PBCs) were applied in all three directions. Cu, Al, Cu-Zr, Si, and ppl-Si liquids each contains 128 000 atoms; the system sizes are sufficiently large to minimize the artifacts associated with the PBC. Two different ensembles (NPT and NVT) were applied to prepare liquid configurations at each temperature of interest. The liquid samples were equilibrated at near or above the melting temperature, under Nose-Hoover thermostat (for NPT and NVT) and barostat (for NPT) with zero external pressure.28 The density for each sample under NVT ensemble was kept constant (the zero pressure density at the melting temperature). For each ensemble, 50 independent instantaneous liquid configurations (ILCs) were recorded and averaged. The corresponding inherent structures (ISs) were obtained by energy minimization using the conjugate gradient method. The GST liquid configurations were obtained from ab initio molecular dynamics simulations using the Vienna ab initio simulation package (VASP),29 with canonical (NVT) ensembles and temperature controlled by the Nose thermostat. Projector augmented-wave (PAW) method with the Perdew-Wang exchange-correlation potentials were adopted.34, 35 The simulation was performed at the point only. The GST samples each had 378 atoms (84 Ge atoms, 84 Sb atoms, and 210 Te atoms under PBCs) and were melted at 3000 K for 2000 time steps (6 ps) and then quenched to a specific temperature of interest, followed by relaxation. Six independent configurations were then selected for structural analysis. To investigate the distance-dependent dilation strain between the instantaneous liquid configurations and the corresponding inherent structure, we define σ thermal (r) following.36–38 By grouping the atomic pairs within the shell between (r − r) and (r + r) in the instantaneous liquid configurations, we solve the affine transformation matrix Jr which best map the following:
Jr rjIiLC → rjIiS ; r − r < rj0i ≤ r + r.
FIG. 1. PDFs of Al liquids at several temperatures under zero external pressure (NPT ensemble). The inset shows a magnified view of their first peaks.
yy
ηxx +η +ηzz
σthermal (r) = r 3r r . In our definition, positive σ thermal (r) means that the shell at the distance of r in the ILCs will dilate upon energy minimization to reach the corresponding IS. III. RESULTS AND DISCUSSIONS
First, we examine the temperature dependence of PDFs in metallic melts. In Fig. 1(a), PDFs of ILCs for Al melts from 900 K to 1900 K under zero external pressure (NPT ensemble) are plotted. The first peak, blown-up in the inset, exhibits the “anomalous” behavior noticed in a number of previous experiments and simulations:12–16 the position corresponding to the maximum peak height shifts to shorter distance with increasing temperature. We also simulated liquids with a fixed density. Figure 2 presents the PDFs of both ILC and IS of Al melts under the constant-volume (NVT ensemble) simulation. As expected, the PDFs of the ISs at various temperatures all collapse onto one single profile, exhibiting
(3)
The Lagrangian strain matrix can then be calculated as ηr =
1 Jr JrT − I . 2
(4)
FIG. 2. PDFs of Al liquids and their inherent structures at various temperatures (simulated under NVT ensemble, i.e., the density is fixed and temperature independent).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 66.173.107.86 On: Wed, 12 Mar 2014 20:35:03
064501-3
Ding et al.
FIG. 3. Schematic showing instantaneous liquid configurations (solid circles) and the corresponding inherent structures (open circles), from the perspective of potential energy landscape (PEL).
negligible temperature dependence. This is consistent with the PEL theory, i.e., equilibrium liquids at a constant density are iso-configurational.22, 30, 31 With no volume expansion nor evolution of inherent structure, here the first peak position in PDFs of ILCs clearly shifts to the left with increasing temperature (similar to the NPT results in Fig. 1 and Ref. 22). As such, Fig. 2 demonstrates that the temperature-dependent shifts in the PDFs of the ILC can arise from thermal effects alone, in the absence of any appreciable evolution of the IS22, 30 or positive/negative thermal expansion.
J. Chem. Phys. 140, 064501 (2014)
To illustrate the difference between the ILC and its corresponding IS state, Fig. 3 schematically depicts a PEL,30, 31 with solid and open circles representing the ILCs and ISs, respectively. The IS is the local minimum of the potential energy basin in which the ILC resides, after removing the temperature effects through energy minimization (the IS may be regarded as the go-to configuration at zero temperature). From this perspective, the results in Fig. 1 indicate that it is the increased thermal excitation at higher temperatures that induces a re-distribution of atoms over a wider range of distances to cause an anomalous shift of the first peak maximum towards the shorter r side. This latter conjecture can be shown to hold, in the absence of changes in system IS/density, as done in Fig. 2. Not all liquids behave this way, though. In Fig. 4, the PDFs of ILC and the corresponding IS for (a) Cu, (b) Cu64 Zr36 , (c) Si, and (d) Ge2 Sb2 Te5 39 are presented. For Cu and Cu64 Zr36 ILCs, their first peak positions in PDFs shift to shorter r with increasing temperature, just like Al (see Figs. 1 and 2) and other metallic liquids. On the contrary, the first peak position in PDFs of Si and Ge2 Sb2 Te5 liquids moves to larger r with rising temperature, which agrees with other experimental and simulation data of molten covalent Si, Ge, and polymers.22–26 Meanwhile, the PDFs of IS for each of the liquid systems in Fig. 4 always overlap (there are
FIG. 4. PDFs of (a) Cu, (b) Cu64 Zr36 , (c) Si, and (d) GST liquids and their corresponding inherent structures at various temperatures at constant density (NVT ensemble).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 66.173.107.86 On: Wed, 12 Mar 2014 20:35:03
064501-4
Ding et al.
minor variations for GST liquids due to the limited box-size and number of configurations in DFT calculation). Also, by comparing the PDFs of ILC and IS in Fig. 4 (and Fig. 2), it can be seen that for metallic liquids, the first peak position in PDFs of ILC is on the left-hand (lower r) side relative to the peak position in the IS PDFs (note that the splitting of the first peak in Cu-Zr inherent structure comes from the different atomic radius of Cu and Zr atoms), but for Si and GST liquids the peak position of the ILC PDF is on the right-hand (larger r) side. Since IS can be regarded as a configuration at zero temperature, the above is also consistent with the observation that with increasing temperature, the ILC peak position shifts towards the left for metallic liquids, but to the right for covalent liquids. In other words, the position corresponding to the IS corresponds to the trend extrapolated to zero temperature. In liquids, the nearest neighbors (NNs) and CNs are not unambiguously defined. In order to evaluate NN or CN, a distance cutoff at the minimum after the first peak in PDF is widely employed. But caution should be exercised when it is used to analyze subtle changes in NN or CN. The distance cutoff is not set at zero PDF intensity, and the minimum (the dip) in fact corresponds to a fairly high intensity (see Fig. 1), especially when considering that for high-temperature liquids the peak intensity distribution spans a rather wide r range. In our constant-density simulations, since ILC at all temperatures is iso-configurational with virtually identical IS, the CNs of their IS should also be the same no matter which method is employed to determine the NNs. In contrast, if we analyze ILC using the distance cutoff at first minimum, the CN will decrease with increasing temperature (not shown here), as observed in Ref. 12. This is because from IS to ILC, and from ILC at a lower temperature to that at a higher temperature, the profile (shape) of the PDF peak changes (broadens) due to thermal activation, such that the pre-set cutoff misses part of the first neighbors that have moved to longer r. Therefore, the apparent change in CN for the ILCs, when determined by the distance cutoff method, is merely a dynamical effect due to thermally induced spreading. The same is true for the “first peak position,” as will be explained next. In Fig. 4, it is observed that the profile of the first peak in PDFs becomes increasingly asymmetric with increasing temperature. To evaluate the asymmetry, we employ the ratio of rB /rA , with the two distances rB and rA defined at the full width of half maximum height of the first PDF peak, see inset in Fig. 5 for a Cu liquid at 1300 K. For a symmetric peak, rB /rA will be equal to 1. Fig. 5 shows that the rB /rA for Al, Cu, and Cu-Zr liquids (Si and GST liquids are excluded for this analysis because their first peaks are not sufficiently sharp) is always greater than 1 and increases with increasing temperature. This clearly indicates that higher temperature leads to higher asymmetry of the first peak in PDFs for equilibrium metallic liquids. Specifically, for the first PDF peak the curve at shorter distances (the left wing of the peak) is steeper than that at larger distances (right wing). This is due to the shape of the typical interatomic interaction potential for pair interactions (the potential energy versus r curve):20, 22 the anharmonicity of the potential grows as atoms reach higher energy
J. Chem. Phys. 140, 064501 (2014)
FIG. 5. Temperature dependence of the asymmetry of the first peak in PDFs (evaluated by rB /rA ) for Cu, Al, and Cu-Zr liquids (constant density). The inset schematically shows the definition of rB and rA (at the full width of half maximum of the first peak for Cu liquid at 1300 K; purple circles).
levels at elevated temperatures. More discussions on this point will be presented later. To see the re-distribution due to temperature, we compare the ILC and the IS, at a constant density, to observe how much the atoms move upon energy minimization. The σ thermal (r) defined in Sec. II can measure the distance-dependent dilation strain from ILC to the corresponding IS. For instance, a positive σ thermal (r) means that the shell at the distance of r in ILC will dilate when reaching the corresponding IS and vice versa. As presented in Fig. 6(a), σ thermal (r) of Cu liquids at 1300 K and 2500 K are very similar and the dilation strain ranges between −6% and 21%. In other words, some interatomic distances would dilate as much as 21% upon energy minimization, while some others shrink by a few percent. Also, a zero value of σ thermal (r) in Fig. 6(a) means that a balance is reached between dilation and contraction for that shell at the distance of r. For comparison, Fig. 6(b) shows the average r number of atomic pairs within the distance of r, N (r) = 0 4π r 2 ρ0 g(r)dr(where ρ 0 is the average number density) for Cu liquids at 1300 K, 2500 K, and IS (see left axis), as well as the PDF of IS for Cu liquid (see right axis). The dashed lines in Fig. 6 indicate the locations where σ thermal (r) is equal to zero. We observe that (i) σ thermal (r) = 0 corresponds to the crossing between the N(r) of IS and the N(r) at high temperatures and (ii) the zero value of σ thermal (r) also corresponds to either the peak or valley position (except R2 ) in the PDF of (r) < 0 corresponds IS. Specifically, σ thermal (r) = 0 and dσthermal dr (r) >0 to peak position, whereas σ thermal (r) = 0 and dσthermal dr coincides with the minimum valley in PDF in Fig. 6(b). We observe that the position of R2 does not directly point to peak position; this is due to the splitting of this peak. But R2 /R1 ≈ 1.8 is consistent with the prediction by spherical-periodic ,40 where n is an integer, order (SPO), RnSP O /R1SP O = 15 + 4n 5 as well as the average position of the second atomic shell as calculated for five different amorphous alloys.41 Therefore, R2 , as an average position of the second shell, is in line with the interpretation of R1 and R3 (R3 /R1 ≈ 2.6, also in agreement with SPO prediction). In sum, the observations above
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 66.173.107.86 On: Wed, 12 Mar 2014 20:35:03
064501-5
Ding et al.
FIG. 6. (a) Distance dependence of the dilation strain (σ thermal ) from instantaneous liquid configuration to the corresponding inherent structure in Cu liquids at 1300 K and 2500 K. The dashed lines indicate the position r with zero value of σ thermal . (b) N (r) = 0 4π r 2 ρ0 g(r)dr for Cu at 1300 K, 2500 K and the corresponding inherent structure (left axis), and the PDF of inherent structure of Cu liquids at 1300 K (right axis).
provide a general picture of the mapping between an ILC and its corresponding IS: (i) in ILC, the interatomic distances in IS spread out to become both shorter and longer under thermal agitation; (ii) the peaks (centers of atomic shells) and minima (valleys separating two neighboring shells) reach the specific condition of σ thermal (r) = 0 under thermal agitation. From the information above, two reasons come to mind that can explain why the position of the maximum intensity of the first peak in PDFs of metallic melts shifts to shorter r with increasing temperature. First, as shown in the preceding paragraph and Fig. 6, thermal agitation of atoms in metallic liquids leads to both shorter and longer bonds (i.e., negative and positive departure from the bond length in IS, RIS , which corresponds to the IS first peak position). The average bond length (Rave ), estimated considering CN to be the same as that in the corresponding IS, is plotted in Fig. 7(a) as a function of temperature for each of the two groups of atoms, one with shorter bonds (r < RIS ) and the other with longer ones (r > RIS ). Starting from the two data points at 0 K for IS, we observe that with increasing temperature, the two groups diverge, indicating an expanding atomic distance range. The Rave of the r > RIS group is larger than that of the r < RIS group. This, again, is because typical pair interactions are asymmetric and anharmonic, which becomes especially pronounced at high temperatures. The potential-distance curve is usually steeper on the shorter distance side but becomes much less so for longer r. In other words, the resistance the thermally excited atoms have to combat with is stiffer on the r < RIS side but lower on the r > RIS side. As a result, under thermal activation it is easier for atoms to move further away from the
J. Chem. Phys. 140, 064501 (2014)
central atom to reach r > RIS , with a distance range wider than that on the r < RIS side. To counter balance these latter atoms such that the system average is still at RIS , a higher density of atoms is expected to populate on the narrower r < RIS side. This asymmetric atomic distribution is illustrated in Fig. 7(b), reflected by a higher PDF peak intensity at smaller r in the left wing, together with broadening towards larger r on the right side. One can see this by hypothetically assuming the same number of atoms drifting to the left of RIS and to the right of RIS . Those on the right tend to take r positions further away from RIS , while those on the left reside closer to RIS over a smaller r range. Consequently the number density appears to be higher on the smaller r side. This is the reason for the enhanced asymmetry of the first PDF peak at high temperatures (also seen in Fig. 5). Such a skewing of the peak results in an apparent shift of the peak maximum position to lower r with increasing temperature, even when the system has the same average bond length (RIS ) in our constant-volume (and constant IS and CN) simulation. A similar trend is expected under constant pressure conditions. Second, another contributor to the change of the first peak position/shape is the normalization of 4π r2 in the calculation of g(r). As g(r) = 4πrS(r) 2 ρ dr , where S(r) is the average 0 number of atoms within the shell with atomic distance between (r − 12 dr) and (r + 12 dr), the g(r) peak position will be different from that in the radial distribution function, RDF (= 4π r2 g(r)), as seen from the difference between RPDF and RRDF in Fig. 7(c). This difference between g(r) and 4π r2 g(r) increases with rising temperature (see the plot of RRDF -RPDF in Fig. 7(d)). For a given change in S(r) due to thermal motion, if it happens on the longer r side it would drop the PDF intensity more, when compared with the effect of the same change in S(r) on shorter r side. This is because RDF at higher r is divided by a larger normalization factor 4π r2 . This causes a relative enhancement of the intensity of the peak on the leftwing side (smaller r), contributing to the peak asymmetry and the apparent position shift of the maximum intensity of the peak. The last question to address is the origin of the different behavior observed in melts of covalent Si and GST. In Fig. 8, we plot the distance-dependent N(r) for Cu and Si at different temperatures. The average first peak positions (RIS ) in PDFs are chosen as the reference point. A clear difference between Cu and Si can be observed at the position of RIS and in the region immediately below RIS : for Cu we see more atoms than IS, whereas for Si there are fewer atoms than IS. In other words, there appears to be an additional driving force in liquid Si to move a higher fraction of surrounding nearest neighbor atoms to distances a bit longer than RIS than in liquid Cu. This is the reason why the position of the maximum of the first peak in the PDFs of Si (and GST and other covalent liquids) shifts to longer r with increasing temperature, which contrasts with the case in metallic melts such as Cu, Al, and Cu-Zr. This difference is rooted in an extra contribution from the directional covalent bonding. To see this we compare the standard Stillinger-Weber potential22 and a modified version without the three-body (directional) term, as described in Sec. II. The potential with pair interaction only is plotted in
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 66.173.107.86 On: Wed, 12 Mar 2014 20:35:03
064501-6
Ding et al.
J. Chem. Phys. 140, 064501 (2014)
FIG. 7. (a) Temperature dependence of the average interatomic distance for atom pairs with r ≤ RIS (shorter distance) and r ≥ RIS (longer distance) in Cu melts; (b) schematic showing increasing asymmetry of the first PDF peak with increasing temperature; (c) pair distribution function (PDF) compared with radial distribution function (RDF) for Cu liquid at 1300 K, as well as the corresponding ISs; and (d) temperature dependence of the difference between the first peak positions in RDF and PDF (denoted as RRDF − RPDF ).
the inset of Fig. 9 (solid line). This pair-potential simulated liquid without the three-body angular term is called ppl-Si hereafter for convenience, and the temperature dependence of the PDF of ppl-Si in reduced units is shown in Fig. 9. In contrast with that in Si, the first peak now also shifts to shorter r with increasing temperature, resembling the behavior in metallic liquids. This result clearly demonstrates the
role of the directionality of the covalent bonds. Specifically, liquid Si originates from covalent Si. In the melt the bonding becomes increasingly more metallic, while retaining certain degree of covalency. The liquid structure exhibits a more pronounced signature of ∼60◦ bond angle, with increasing deviation from its equilibrium bond angle (corresponding to cos (θ ) = 1/3).42 As seen in Eq. (2), the interaction energy would be
FIG. 8. Distance-dependent N(r) for ILC and IS of (a) Cu and (b) Si liquids. Colored regions are for r ≤ RIS . Note that when the number of neighbors is counted up to RIS (corresponding to N(r) at r = RIS ), the number of neighbors in the ILC is larger than that in the IS for Cu, but the opposite is true for Si.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 66.173.107.86 On: Wed, 12 Mar 2014 20:35:03
064501-7
Ding et al.
FIG. 9. PDFs of ppl-Si (potential modified to eliminate the three-body angular term) at various temperatures (in reduced units). The inset compares the modified potential (solid line) with the original Stillinger-Weber potential of Si.
elevated by the angular term even for unchanged bond length. For a simplistic configuration of only three atoms (see Fig. 9), which form an equilateral triangle with bond length r, their reduced potential-distance curve is plotted as dotted line in the inset of Fig. 9 if the angular term in Eq. (2) is included. A comparison with the curve from pair potential only (solid line) indicates that the angular term in Stillinger-Weber potential results in a shift of the potential-distance curve to higher r and makes it shallower on the high r side when compared with the pair interaction case. This moves some nearest neighbors to slightly longer distances. In other words, thermal effect in Si liquids enhances the deviation of the bond angle, such that atomic bonds tend to stretch to longer lengths. As a result of these atoms accumulated at such distances, the peak maximum in PDFs shifts to the right (longer r) with increasing temperature, in contrasts with the behavior in metallic liquids.
IV. CONCLUSION
In summary, we have employed classical and ab initio molecular dynamics simulations to systematically study the effects of temperature/thermal excitation on the atomic pair distributions in several representative liquids (both metallic and covalent). A comparative analysis of the ILC and IS has been conducted. Our results show that: (1) The seemingly anomalous peak shift with increasing temperature observed for metallic liquids is not truly anomalous, but intrinsic, reproducible, and predictable, see (2) below. It can happen even when the liquid experiences no change in inherent topological12 or chemical (such as in alloys with two or more constituent elements) short-range ordering. (2) One should study the entire peak profile, rather than just the position corresponding to the maximum intensity. The latter itself is not a meaningful and reliable metric to monitor and understand the changes in CN, bond length (interatomic distance), and local order. The
J. Chem. Phys. 140, 064501 (2014)
PDF as measured includes both configurational (inherent structure, local minimum in PEL) and thermal (deviation from the PEL minimum due to thermal activation) contributions. One intrinsic reason for the position/shape evolution is the thermally induced asymmetric redistribution of neighboring atoms to shorter and longer distances due to the anharmonic nature of the interatomic interaction potential. With increasing temperature, the PDF peak (the distribution of the atoms in the first shell) becomes increasingly more asymmetric. It would be inadvisable to draw conclusions based solely on a comparison of the positions of the intensity maximum of two peaks with different degree of asymmetry. (3) The opposite trend of first peak shift for metallic and covalent systems is systematic (although exceptions may exist, such as multi-component liquids or supercooled liquids under structural ordering). In particular, we note that the directionality term in the interaction potential influences the population and spreading of atomic distributions under thermal agitation, hence the direction and degree of peak shift with increasing temperature.
ACKNOWLEDGMENTS
J.D. and E.M. were supported at JHU by U.S.-DOEBES, Division of Materials Sciences and Engineering, under Contract No. DE-FG02-09ER46056. The computer simulations were performed using the NERSC supercomputers. Y.Q.C. was supported by the Scientific User Facilities Division, Office of Basic Energy Sciences, U.S. Department of Energy. 1 P.
G. Debenedetti and F. H. Stillinger, Nature (London) 410, 259 (2001). C. Frank, Proc. R. Soc. London, Ser. A 215, 43 (1952). 3 C. A. Angell, MRS Bull. 33, 544–555 (2008). 4 Z. H. Jin, P. Gumbsch, K. Lu, and E. Ma, Phys. Rev. Lett. 87, 055703 (2001). 5 Y. T. Shen, T. H. Kim, A. K. Gangopadhy, and K. F. Kelton, Phys. Rev. Lett. 102, 057801 (2009). 6 D. Holland-Moritz, D. M. Herlach, and K. Urban, Phys. Rev. Lett. 71, 1196 (1993). 7 J. C. Bendert, A. K. Gangopadhyay, N. A. Mauro, and K. F. Kelton, Phys. Rev. Lett. 109, 185901 (2012). 8 C. Tang and P. Harrowell, Nat. Mater. 12, 507–511 (2013). 9 J. Ding, Y. Q. Cheng, H. W. Sheng, and E. Ma, Phys. Rev. B 85, 060201 (2012). 10 J. Ding, Y. Q. Cheng, and E. Ma, Acta Mater. 61, 3130–3140 (2013). 11 Y. Q. Cheng and E. Ma, Prog. Mater. Sci. 56, 379–473 (2011). 12 H. B. Lou, X. D. Wang, Q. P. Cao, D. X. Zhang, J. Zhang, T. D. Hu, H. K. Mao, and J. Z. Jiang, Proc. Natl. Acad. Sci. U.S.A. 110, 10068–10072 (2013). 13 G. W. Lee, A. K. Gangopadhyay, K. F. Kelton, R. W. Hyers, T. J. Rathz, J. R. Rogers, and D. S. Robinson, Phys. Rev. Lett. 93, 037802 (2004). 14 P. Ganesh and M. Widom, Phys. Rev. B 74, 134205 (2006). 15 T. H. Kim and K. F. Kelton, J. Chem. Phys. 126, 054513 (2007). 16 N. A. Mauro, J. C. Bendert, A. J. Vogt, J. M. Gewin, and K. F. Kelton, J. Chem. Phys. 135, 044502 (2011). 17 N. A. Mauro, V. Wessels, J. C. Bendert, S. Klein, A. K. Gangopadhyay, M. J. Kramer, S. G. Hao, G. E. Rustan, A. Kreyssig, A. I. Goldman, and K. F. Kelton, Phys. Rev. B 83, 184109 (2011). 18 N. A. Mauro, W. Fu, J. C. Bendert, Y. Q. Cheng, E. Ma, and K. F. Kelton, J. Chem. Phys. 137, 044501 (2012). 19 K. Georgarakis, D. V. Louzguine-Luzgin, J. Antonowicz, G. Vaughan, A. R. Yavari, T. Egami, and A. Inoue, Acta Mater. 59, 708–716 (2011). 20 W. Kob and H. Andersen, Phys. Rev. E 51, 4626 (1995). 2 F.
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 66.173.107.86 On: Wed, 12 Mar 2014 20:35:03
064501-8 21 T.
Ding et al.
H. Kim, G. W. Lee, B. Sieve, A. K. Gangopadhyay, R. W. Hyers, T. J. Rathz, J. R. Rogers, D. S. Robinson, K. F. Kelton, and A. I. Goldman, Phys. Rev. Lett. 95, 085501 (2005). 22 F. H. Stillinger and T. A. Weber, Phys. Rev. B 31, 5262 (1985). 23 S. Sastry and C. A. Angell, Nat. Mater. 2, 739–743 (2003). 24 A. Filipponi and A. D. Cicco, Phys. Rev. B 51, 12322 (1995). 25 K. S. Chang, Y. L. Wang, C. H. Kang, H. J. Wei, Y. H. Weng, and K. L. Tung, J. Membr. Sci. 382, 30 (2011). 26 O. Alexiadis and V. G. Mavrantzas, Macromolecules 46, 2450–2467 (2013). 27 D. C. Rapaport, The Art of Molecular Dynamics Simulation (Cambridge University Press, 2004). 28 M. P. Allen and D. J. Tidesley, Computer Simulation of Liquids (Clarendon Press, Oxford, 1989). 29 G. Kresse and J. Furthmüller, J. Comput. Mater. Sci. 6, 15–50 (1996). 30 F. H. Stillinger, Science 267, 1935–1939 (1995).
J. Chem. Phys. 140, 064501 (2014) 31 D.
Wales, Energy Landscapes: Applications to Clusters, Biomolecules and Glasses, 1st ed. (Cambridge University Press, 2004). 32 H. W. Sheng, M. J. Kramer, A. Cadien, T. Fujita, and M. W. Chen, Phys. Rev. B 83, 134118 (2011). 33 Y. Q. Cheng, E. Ma, and H. W. Sheng, Phys. Rev. Lett. 102, 245501 (2009). 34 P. E. Blöchl, Phys. Rev. B 50, 17953–17979 (1994). 35 G. Kresse and D. Joubert, Phys. Rev. B 59, 1758–1775 (1999). 36 M. L. Falk and J. S. Langer, Phys. Rev. E 57, 7192 (1998). 37 J. Li, Modell. Simul. Mater. Sci. Eng. 11, 173 (2003). 38 J. Ding, Y. Q. Cheng, and E. Ma, Appl. Phys. Lett. 101, 121917 (2012). 39 M. Xu, Y. Q. Cheng, H. W. Sheng, and E. Ma, Phys. Rev. Lett. 103, 195502 (2009). 40 P. Häussler, Phys. Rep. 222, 65 (1992). 41 J. Ding, Y. Q. Cheng, and E. Ma, “Atomic order beyond short-range scale in metallic glasses and liquids” (unpublished). 42 C. S. Liu, Z. G. Zhu, J. C. Xia, and D. Y. Sun, Phys. Rev. B 60, 3194 (1999).
This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 66.173.107.86 On: Wed, 12 Mar 2014 20:35:03
