Lattice Dynamics and Thermodynamics of Molybdenum from First-Principles Calculations
of the ultrahigh pressure scale. The equation of state (EOS) of Mo at high pressure is being used as a calibration standard to the ruby fluorescence in diamond anvil compression (DAC) experiments.1 Because of its important position in the field of material science and condensed matter science, Mo has attracted tremendous experimental and theoretical interest in its wide range of properties recently. At ambient condition, Mo is in body-centered-cubic (bcc) structure and melts at 2890 K.2 But what is the most stable phase of Mo under ultrahigh pressure? Experimentally, the shock wave (SW) acoustic velocity measurements showed that there was a sharp break on the Hugoniot curve at about 210 GPa (at a calculated temperature of 4100 K), which indicated that a solid-solid phase transition occurred prior to melting at 390 GPa (at a calculated temperature of 10 000 K).3 To compare with the SW data, Vohra and Ruoff investigated the static compression of Mo by energy-dispersive X-ray diffraction and found that the bcc Mo was stable up to 272 GPa at 300 K.4 The phase transition in shock compression at 210 GPa was not observed. By further X-ray diffraction investigation, Ruoff et al. showed that Mo was stable in a bcc structure up to at least 560 GPa at room temperature.5 Theoretically, Moriarty suggested that Mo was stable in the bcc structure up to 420 GPa, where it transformed to a hexagonal close-packed (hcp) structure and then at 620 GPa to a face-centered closepacked (fcc) structure.6 Later, Boettger7 and Christensen et al.8 showed that the hcp phase of Mo was not stable, and the bcc phase transformed directly to the fcc phase at 700 GPa.
Belonoshko et al. confirmed the results of Boettger and showed the transition pressure was 720 GPa at zero pressure.9 By calculating the Gibbs free energies of the bcc and fcc Mo in the pressure range from 350 to 850 GPa at room temperatures up to 7500 K, Belonoshko et al. found that Mo had lower free energy in the fcc phase than in the bcc phase at elevated temperatures.10 Shortly after these new results were reported, Cazorla et al. found that the hcp Mo was noticeably more stable above 350 GPa at high temperature by calculating the Helmholtz and Gibbs free energies of the bcc, fcc, and hcp phases.11 The other intriguing problem is the melting properties. For the transition metals, such as Mo, Ta, and W, there are enormous discrepancies in melting curves between laser-heated DAC12-17 and SW3,18 methods. As for Mo (as well as Ta and W), several thousand degrees of discrepancies exist in extrapolating from DAC pressures of around 100 GPa12-15,17 to SW pressure of 390 GPa.18 As is known that the overestimation of the melting temperature exists in SW experiments, Errandonea13 corrected the SW data by considering 30% superheating. The revised melting temperatures are located at 7700 ( 1500 K (390 GPa), also much larger than the melting temperature (just above 4000 K) at this pressure extrapolated from DAC experiments. Results from the empirical and phenomenological melting models are dependent on the selection of the parameters.13,19-22 It is shown that the choice of different sets of parameters leads to huge differences in the melting temperatures at high pressure. In addition, theoretical results are consistent with the SW data at high pressure but diverge from DAC data below 100 GPa. All these results are inadequate to explain the extreme discrepancies in extrapolating the DAC data to the SW data.9-11,23 The highpressure melting curve of Mo still remains inconclusive up to now.
The low-temperature phonon spectrum of Mo as measured in inelastic neutron scattering experiments exhibited a variety of anomalies: the large softening near the H point, the T2 branch at the N point, and the longitudinal branch did
not display the rounded dip near q ) 2/3 [111], which was typical for "regular" monatomic bcc metals.24-26 With temperature increasing, the H point phonon displayed anomalous stiffening, which had been proposed to arise from either intrinsic anharmonicity of the interatomic potential or electron-phonon coupling.26 Theoretically, using the molecular dynamics (MD) simulations with environment-dependent tight-binding parametrization, Haas et al. reproduced the weakening of the phonon anomalies as the temperature increased.27 Later, Farber et al. determined the lattice dynamics of Mo at high pressure to 37 GPa using high-resolution inelastic X-ray scattering (IXS).28 Meanwhile, they calculated the quasiharmonic phonon spectrum up to the highest experimental pressure by linear response theory. Both the experimental and theoretical results showed an obvious decrease in the relative magnitude of the H point phonon anomaly under compression. Recently, Cazorla et al. obtained the phonon dispersion curves of the bcc Mo at zero pressure using the small displacement method,29 but under larger compression, there are no experimental and theoretical studies. The first-principles density functional theory is very successful in predicting the high-pressure behavior of phonon dispersion relations and their concomitant anharmonic effects.30 It is more effective to connect the phonon properties directly to the lattice dynamics under pressure, temperature, or their combination. One of the main purposes of this work is to investigate the lattice dynamics and thermodynamics of Mo under high pressure and temperature.
The aim of the present work is multiple. First, density functional theory with the generalized gradient approximation (GGA) has been used for first-principles studies of the phase transition and elastic properties of Mo under high pressure, and then the quasiharmonic approximation (QHA) has been applied to the study of the lattice dynamic properties, the thermal EOS, and the thermodynamic properties. The organization of this paper is as follows. In Section 2, we give a brief description of the theoretical computational methods. The results and detailed discussions are presented in Section 3. A short conclusion is drawn in the last section.
Computational Methodology
For many metals and alloys, the Helmholtz free energy F can be accurately separated as
F(V, T) = Estatic(V) + Fphon(V, T) + Felec(V, T) (1)
where Estatic(V) is the energy of a static lattice at zero temperature; Felec(V,T) is the thermal free energy arising from electronic excitations; and Fphon(V,T) is the phonon contribution. Both Estatic(V) and Felec(V,T) can be obtained from first-principles calculations directly. Density functional perturbation theory (DFPT) is a well-established method for calculating the vibrational properties from first-principles in the framework of QHA.31,32 Including part of the anharmonic effects by considering the volume dependence of phonon frequencies gives access to the thermal expansion, thermal EOS, and thermodynamic properties. Our calculations were performed within the GGA, as implemented in the QUANTUM-ESPRESSO package.33 A nonlinear core correction to the exchange-correlation energy function was introduced to generate Vanderbilt ultrasoft pseudopotential for Mo with the valence electrons configuration of 4d55s1. In addition, the pseudopotential was generated with a scalar-relativistic calculation using GGA according to the recipe of Perdew-Wang 91.34
During our calculations, we made careful tests on k and q grids, the kinetic energy cutoff, and many other parameters to guarantee phonon frequencies and free energies to be well converged. Dynamical matrices were computed at 29 wave (q) vectors using an 8 × 8 × 8 q grid in the irreducible wedge of the Brillouin zone. The kinetic energy cutoff Ecutoff was 60 Ry,
and the k grids used in total energy and phonon calculations were 20 × 20 × 20 and 14 × 14 × 14 Monkhorst-Pack (MP) meshes.35 The self-consistent calculation was terminated when the total energy difference in two successive loops was less than 10-12 Ry. A Fermi-Dirac smearing width of 0.02 Ry was applied for Brillouin zone integrations in phonon frequency calculations, and in the calculations of static energy and thermal electronic excitations we treated the smearing width of Fermi-Dirac as the physical temperature of electrons. The geometric mean phonon frequency is defined by
where ωqj is the phonon frequency of branch j at wave vector q, and Nqj is the number of branches times the total number of q points in the sum. With the tested parameters, the geometric mean phonon frequency was converged to 1 cm-1. The elastic constants are defined by means of a Taylor expansion of the total energy, E (V, δ), for the system with respect to a small strain δ of the lattice primitive cell volume V. The energy of a strained system is expressed as follows
where E(V0, 0) is the energy of the unstrained system with equilibrium volume V0; τi is an element in the stress tensor; and i is a factor to take care of the Voigt index.36 The independent elastic constants in cubic structure are C11, C12, and C44. For the calculations of three elastic constants, we considered three independent volume-nonconserving strains
The specific energy of the crystal deformed in accord with the matrix kε 1 was calculated as a function of the strain magnitude γ and was used in an equation that can be written as
Results and Discussion
3.1. Phase Transition. We obtained the static energy-volume (E-V) curves for the bcc, fcc, and hcp structures. From the static energy differences shown in Figure 1(a), we can see that the hcp structure is unstable in all ranges of volumes. The bcc structure is stable above 8 Å3. It means that the bcc Mo transforms to the fcc Mo under compression. By fitting the E-V data to the fourth-order finite strain EOS,39,40 we can obtain the enthalpy as a function of pressure. The enthalpy differences are shown in Figure 1(b). We find that Mo is stable in the bcc structure up to 703 ( 19 GPa, and then it transforms to the fcc phase. The hcp structure can not be the stable phase in the whole range of pressures. The equilibrium volume V0, bulk modulus B0, and its pressure derivatives B′ and B′′ are listed in Table 1. The agreement of our results with experimental data4,41 is very good. From Table 1, we can see that the bulk modulus of fcc
Mo at zero pressure is smaller than that of bcc Mo. It seems incompatible because the fcc phase is more densely packed than the bcc phase. Actually, the fcc phase of Mo is not observed at zero pressure and temperature. Comparing the first and second pressure derivatives of bulk modulus of bcc Mo with that of the fcc Mo, it is obvious that the bulk modulus of fcc Mo is larger than that of bcc Mo under ultrahigh pressure. Whether the bcc Mo is stable or not under high temperature and high pressure remains disputable. By calculating the Gibbs free energies of the bcc and fcc Mo from 350 to 850 GPa at room temperature up to 7500 K, Belonoshko et al. found that the bcc Mo transformed to the fcc phase below 700 GPa at high temperature.10 We also calculated the free energy within QHA and obtained the phase diagram of Mo at high pressure and temperature, which is shown as Figure 1(c). It can be seen that our bcc-fcc boundary above 700 GPa is consistent with that obtained by Belonoshko et al. Below 700 GPa, the bcc-fcc boundary from Belonoshko et al. only exists at the pressure range 350-700 GPa. However, from our calculations, the boundary can extend to 232 GPa (4000 K). It is very close to the SW data 210 GPa (4100 K), where the discontinuity in longitudinal sound velocity was detected.3 It may imply that the bcc Mo transforms to the fcc phase at this shock compression. If we extrapolate the bcc-fcc boundary to 119 GPa, the temperature will locate at 3557 K, a little larger than the latest double-sided laser-heated DAC melting temperature 3302 ( 140 K17 at this pressure. So we consider that the DAC melting curve might not be a simple solid-liquid phase boundary but the indication of an unusual transition. Indeed Santamarı´a-Pe´rez et al. also suggested that the melted Mo in DAC experiments was a complex structured liquid (may be tetrahedral and icosahedral structures with short-range order), more akin to the supercooled liquid phase than a simple liquid.17 Cazorla et al. also calculated the free energy of the bcc, fcc, and hcp structures and found that the hcp was the most stable phase at this high pressure and temperature region.11 However, their bcc-hcp boundary is very close to the present (and Belonoshko et al.) bcc-fcc boundary.
3.2. Elastic Properties. For each volume of the unit cell, the elastic constants (C11, C12, C44) of the bcc Mo were deduced from a polynomial fit of the strain energy for particular deformations listed in eqs 5-7. Ten symmetric values of γ in the range ( 5% were used to make the strain energy fit at each strain type. For P, the value was calculated from the static equation of state. The obtained elastic constants C11, C12, and C44 at equilibrium structure parameters are 472.3, 160.4, and 106.0 GPa, respectively. The results of these calculations under pressure are presented in Figure 2, which are in good agreement with the energy-dispersive X-ray diffraction data.42 But Duffy
et al. only measured the elastic constants up to 24 GPa, and they showed the pressure derivatives of elastic moduli (linear fit) were 7.3, 3.3, and 0.5.42 From Figure 2, we can see that the pressure dependence of elastic moduli under high pressure can not be described as linear relations. All the elastic moduli can be described as third polynomials perfectly, which are shown
as follows
C11 ) 472.3 + 6.06564P - 5.67 × 10-3P2 +
6.10951 × 10-6P3 (8)
C12 ) 160.4 + 2.71872P - 2.54 × 10-3P2 +
2.58221 × 10-6P3 (9)
C44 ) 106.0 + 2.06781P - 2.13 × 10-4P2 +
2.25069 × 10-7P3 (10)
6.10951 × 10-6P3 (8)
C12 ) 160.4 + 2.71872P - 2.54 × 10-3P2 +
2.58221 × 10-6P3 (9)
C44 ) 106.0 + 2.06781P - 2.13 × 10-4P2 +
2.25069 × 10-7P3 (10)
The theoretical polycrystalline elastic modulus can be determined from the independent elastic constants. The average isotropic shear modulus G and bulk modulus B of polycrystalline can be calculated according to Voigt-Reuss-Hill approximations. 43 Then the isotropically averaged aggregate velocities can be obtained as follows
VP ) [(B + 4/3G)/F]1/2 (11)
VB ) (B/F)1/2 (12)
VS ) (G/F)1/2 (13)
VB ) (B/F)1/2 (12)
VS ) (G/F)1/2 (13)
where VP, VB, and VS are the compressional, bulk, and shear sound velocities, respectively. The bulk and shear modulus B and G at high pressure are shown in Figure 3(a), and the aggregate velocities are shown in Figure 3(b). As pressure increases, the B, G, VP, VB, and VS increase monotonously.
3.3. Phonon Dispersions. Figure 4(a) shows the obtained dispersion curves of bcc Mo at zero pressure along several highsymmetry directions in the Brillouin zone (BZ) for both transverse acoustical (TA) and longitudinal acoustical (LA) branches, together with the experimental dispersion curves.26,28 It can be seen that our results agree excellently with the
experimental data. Along Γ-H, the theory does not predict the
experimental data. Along Γ-H, the theory does not predict the
overbending at about q ) 0.8 in the TA branch. At the H point, the theoretic value is a littl larger (∼5%). Along Γ-P-H, the agreement in shape as well as absolute energies is very good. Along Γ N, the calculation of TA [110]〈001〉 captured the small softening near the zone boundary successfully. Our phonon branches at zero pressure spread up to 277 cm-1 at the N point. Meanwhile, the dispersion curves at 17.2 GPa are plotted in Figure 4(b). The phonon branches at this pressure spread up to 322 cm-1 THz at the N point. One finds the excellent agreement of this dispersion curve with the results from inelastic X-ray scattering (IXS).28 We repeated the phonon calculations for the other 13 different volumes. Some of the dispersion curves are plotted in Figure 5(a), which shows the well-known phenomenon in solids; i.e., the phonon frequencies increase as volume decreases, except the values along the Γ-P-H direction and Γ-N in the lower TA branch. As pressure increases, the softening of TA [110]〈001〉 along Γ-N becomes more and more obvious. Under ultra compression (V ) 7.69 Å3), the frequencies along P-H in TA branches soften to imaginary frequencies, indicating a structural instability. Actually, from our static structure calculations, the bcc phase transforms to the fcc phase at about 703 GPa (about V ) 8.07 Å3). To understand the stability of Mo, we also calculated the phonon dispersions of the fcc Mo in 14 volumes. Some phonon
dispersions are shown in Figure 5(b). It is obvious that the fcc Mo is unstable at zero pressure (V ) 16.14 Å3), and with the decrease of volume, all the frequencies of LA and TA increase. The calculation predicted the stability of the fcc Mo by promoting the frequencies of the phonons along Γ to X and Γ to L symmetry lines from imaginary to real. Though the phonon dispersions indicate the stability of the fcc Mo at 9.68 Å3, the enthalpy at this compression is also insufficient to overcome the potential barrier of the bcc-fcc phase transition. From the stabilities of the bcc and fcc Mo, one notes that the phase transition might not be induced independently by the dynamic instability. However, the phase transition occurs before the TA modes go to zero frequency, and the mode softening behaviors are related to the particular mechanism which is responsible for the phase transition.
Conclusions
We employed the density functional perturbation theory to investigate the phase transition, elastic properties, lattice dynamical properties, thermal equation of state, and thermodynamic properties of the bcc structure Mo. By comparing the static energy-volume and enthalpy-volume curves for the bcc, fcc, and hcp structures, we found that Mo is stable in the bcc structure up to 703 ( 19 GPa and then transforms to the fcc structure. The hcp structure can not be the stable phase in the whole range of pressures. With the QHA, the bcc-fcc boundary under high temperature and pressure was obtained. The present results qualitatively agree with those reported by Belonoshko et al. but need the confirmation by further experiments. Meanwhile, to obtain the more accurate phase diagram, the anharmonic phonon-phonon interactions and electron-phonon interactions must be considered completely. The calculated elastic constants at low pressure agree well with the experimental data. Under high pressure, the elastic constants can not be described as linear relations, but the third-order polynomial can
describe the functions of elastic constants versus pressure satisfactorily. From the elastic constants, the bulk and shear moduli and the sound velocities are obtained successfully. The calculated phonon dispersion curves accord excellently with experiments. Under pressure, we captured a large softening along H-P in the TA branches. When the volume is compressed to 7.69 Å3, the frequencies along H-P in the TA branches soften to imaginary frequencies, indicating a structural instability, while with pressure increasing, the phonon calculations on the fcc Mo predicted the stability by promoting the frequencies along Γ to X and Γ to L symmetry lines from imaginary to real. We also investigated the thermal properties including thermal EOS, thermal pressure, volume thermal expansion, and Hugoniot properties. From thermal expansion coefficients and heat capacity, we found that the QHA is valid only up to about the melting point at zero pressure, but under pressure, the validity of QHA can be extended to much higher temperature. On one hand, the QHA includes part of anharmonicity by allowing phonon frequencies to vary with volume. On the other hand, the pressure can suppress part of the anharmonicity by strengthening the bondings among atoms and lowering the vibration of atoms. The Debye temperatures calculated from the elastic constant agree well with that obtained from the Debye approximation.
Cesar Hernandez
19.502.806
Electronica de Estados Solidos
¿Te gustaría tener Hotmail en tu móvil Movistar? ¡Es gratis!
No hay comentarios:
Publicar un comentario