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)

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)

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

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.