Lattice Dynamics and Thermodynamics of Molybdenum from First-Principles Calculations

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