We present the ab initio thermoelastic properties of body-centered cubic molybdenum under extreme conditions obtained within the quasi-harmonic approximation including both the vibrational and electronic thermal excitation contributions to the free energy. The quasi-harmonic temperature-dependent elastic constants are calculated and compared with existing experiments and with the quasi-static approximation. We find that the quasi-harmonic approximation allows for a much better interpretation of the experimental data, confirming the trend found previously in other metals. Using the Voigt–Reuss–Hill average, we predict the compressional and shear sound velocities of polycrystalline molybdenum as a function of pressure for several temperatures, which might be accessible in experiments.
I. INTRODUCTION
Molybdenum, as a refractory 4d transition metal in the same group of tungsten, finds several applications, pure or in alloys with other metals, for its high melting point, mechanical properties, and corrosion resistance. Its thermodynamic properties have been studied by several authors, both experimentally and theoretically.1–13
Among the thermodynamic properties, the elastic constants (ECs) are key parameters that determine the mechanical stability, the velocity of sound waves, and the stress response to external strains. In this regard, particularly useful for applications is the pressure and temperature dependence of the ECs. In molybdenum, the temperature-dependent bulk modulus and ECs have been measured by the ultrasonic technique at room pressure.14–17 Data are available almost until melting (2898 K).17 Pressure derivatives of the ECs are known at room temperature,18 and compressional and shear sound velocities in polycrystalline molybdenum have been measured up to 120 kbar at room temperature.19 A simultaneous measurement of the density allow us to derive from these data the bulk and shear moduli at high pressures. However, information on pressure-dependent elasticity at high temperatures is still missing in the literature.
Theoretically, the temperature-dependent ECs of molybdenum have been calculated within the quasi-static approximation (QSA) at room pressure,20 while the pressure-dependent ECs have been calculated by Koči et al.21 at zero temperature. The measurements on polycrystalline molybdenum have been modeled by ab initio calculations,19 and the ECs along the Hugoniot together with the corresponding compressional and shear sound velocities have been calculated within the QSA.10 Due to the time-consuming phonon calculations for deformed configurations of metallic systems needed for the ab initio quasi-harmonic (QHA) ECs, no paper has addressed these quantities for molybdenum so far.
In recent years, a workflow for the calculation of the QHA ECs within density functional theory (DFT) has been fully integrated in the thermo_pw software,22,23 continuously improved, and applied to several solids (aluminum, silicon, BAs, copper, silver, gold, palladium, platinum, and tungsten).24–27 In many metals, it has been found that the QHA predicts the temperature dependence of the ECs accurately, much more than the QSA.23,28 The 0 K values instead are similar and might differ from experiment, with differences that depend on the exchange and correlation functional and are usually within 10%.
In this paper, we apply this technique to molybdenum. We present a comparison of the temperature dependence of the QSA and the QHA ECs and use the calculated QHA adiabatic ECs to predict the temperature and pressure dependence of the bulk and shear moduli of polycrystalline molybdenum as well as its compressional and shear sound velocities, providing a theoretical prediction that could be useful in future investigations of the thermoelastic properties at high pressures and temperatures.
II. THEORY AND COMPUTATIONAL PARAMETERS
The temperature- and pressure-dependent thermodynamic properties and ECs are calculated by the open source software thermo_pw, which has been discussed in previous publications.23–27,34,35
The calculations presented in this work are done by using DFT as implemented in the Quantum ESPRESSO (QE) package.36,37 The exchange and correlation functionals are the LDA38 and the generalized gradient approximations PBEsol39 and PBE.40
We employ the projector augmented wave (PAW) method41 and a plane wave basis with pseudopotentials from pslibrary.42 We use Mo.pz-spn-kjpaw_psl.1.0.0.UPF, Mo.pbesol-spn-kjpaw_psl.1.0.0.UPF, and Mo.pbe-spn-kjpaw_psl.1.0.0.UPF for LDA, PBEsol, and PBE, respectively. These pseudopotentials have the 4s, 4p, 4d, and 5s states in the valence, while the other states are frozen in the core and accounted for by the nonlinear core correction.43 The lattice constants of 14 reference geometries from 4.784 to 6.084 a.u. with LDA, from 4.8162 to 6.1162 a.u. with PBEsol, and from 4.8922 to 6.1922 a.u. with PBE with an interval of 0.1 a.u. between geometries have been chosen to calculate the free energies. For the wave function cutoffs, we use 100, 90, and 120 Ry, while for the charge density, we use 400, 360, and 480 Ry, for LDA, PBEsol, and PBE, respectively. The Fermi surface has been dealt with by the smearing approach of Methfessel and Paxton44 with a smearing parameter σ = 0.02 Ry. With this smearing, the Brillouin zone integrals converge with a 40 × 40 × 40 k-point mesh.
For six reference geometries (with i = 1, 4, 7, 11, 12, 13), temperature-dependent ECs are calculated by three strain types that lead to a body-centered cubic, a centered tetragonal, and a rhombohedral strained lattice. Each strain type is sampled by six strains, from ɛ = −0.0125 to ɛ = 0.0125 with step size δɛ = 0.005. A thicker k-point mesh of 45 × 45 × 45 is employed on strained configurations. Each one of the 108 strained configurations requires calculations of phonon frequencies and electronic density of states. Phonon frequencies are calculated by density functional perturbation theory (DFPT)45,46 getting the dynamical matrices on a 8 × 8 × 8 q-point grid. These dynamical matrices have been Fourier-interpolated on a 200 × 200 × 200 q-point mesh to evaluate the free-energy and the thermodynamic quantities. The calculations are all performed on the Leonardo supercomputer at CINECA with a GPU version of thermo_pw that optimizes some routines of QE for problems with dense k-points sampling in metallic systems.47
III. RESULTS AND DISCUSSION
In Table I, the equilibrium lattice constants, bulk moduli, and pressure derivatives of the bulk moduli of molybdenum obtained as parameters of a fourth-order Birch–Murnaghan interpolation of the static energy U(V) are listed together with a few selected values from previous calculations and experiment. Our PAW LDA, PBEsol, and PBE values of the lattice constant differ by less than 0.1% from the all-electron values reported in Ref. 29. With respect to experiment (a = 5.936 a.u. at 0 K), the LDA, PBEsol, and PBE errors are −0.9%, −0.3%, and 0.7%, respectively, with LDA and PBEsol below experiment and PBE above. For the bulk modulus, these errors become 10% (LDA), 7% (PBEsol), and −1% (PBE) with respect to the 0 K value 2653 kbar.14
The equilibrium lattice constants (a0), the bulk moduli (BT), and the pressure derivatives of the bulk moduli of molybdenum calculated in this work compared with previous calculations and with experiment.
. | . | T (K) . | a0 (a.u.) . | BT (kbar) . | . |
---|---|---|---|---|---|
This study | LDA | 0 | 5.884 | 2949 | 4.00 |
295 | 5.894 | 2874 | 4.09 | ||
PBEsol | 0 | 5.916 | 2826 | 4.02 | |
295 | 5.926 | 2753 | 4.10 | ||
PBE | 0 | 5.975 | 2617 | 4.08 | |
295 | 5.986 | 2543 | 4.16 | ||
Calc.9 | PW91 | 0 | 5.988 | 2666 | 4.42 |
Calc.10 | PBE | 0 | 5.996 | 2633 | 4.21 |
Calc.21 | PBE | 0 | 6.001 | 2610 | 4.5 |
Calc.12,a | PBE | 0 | 5.992 | ||
Calc.29 | LDA | 0 | 5.888 | ||
PBEsol | 0 | 5.920 | |||
PBE | 0 | 5.979 | |||
Calc.2 | LDA | 0 | 5.880 | 3010 | 3.99 |
298b | 5.891 | 2950 | 4.01 | ||
PBEsol | 0 | 5.914 | 2870 | 4.02 | |
298b | 5.925 | 2800 | 4.05 | ||
PBE | 0 | 5.981 | 2620 | 4.14 | |
298b | 5.993 | 2550 | 4.17 | ||
Model4 | 300 | 5.945 | 2600 | 4.21 | |
Model30 | 300 | 5.944 | 2605 | 4.05 | |
Expt.2,a | 300 | 5.944 | 2610 | 4.06 | |
Expt.1 | 300 | 5.951 | 2608 | 4.46 | |
Expt.14 | 0 | 2653 | |||
Expt.31 | 2610 | 4.65c/3.95d |
. | . | T (K) . | a0 (a.u.) . | BT (kbar) . | . |
---|---|---|---|---|---|
This study | LDA | 0 | 5.884 | 2949 | 4.00 |
295 | 5.894 | 2874 | 4.09 | ||
PBEsol | 0 | 5.916 | 2826 | 4.02 | |
295 | 5.926 | 2753 | 4.10 | ||
PBE | 0 | 5.975 | 2617 | 4.08 | |
295 | 5.986 | 2543 | 4.16 | ||
Calc.9 | PW91 | 0 | 5.988 | 2666 | 4.42 |
Calc.10 | PBE | 0 | 5.996 | 2633 | 4.21 |
Calc.21 | PBE | 0 | 6.001 | 2610 | 4.5 |
Calc.12,a | PBE | 0 | 5.992 | ||
Calc.29 | LDA | 0 | 5.888 | ||
PBEsol | 0 | 5.920 | |||
PBE | 0 | 5.979 | |||
Calc.2 | LDA | 0 | 5.880 | 3010 | 3.99 |
298b | 5.891 | 2950 | 4.01 | ||
PBEsol | 0 | 5.914 | 2870 | 4.02 | |
298b | 5.925 | 2800 | 4.05 | ||
PBE | 0 | 5.981 | 2620 | 4.14 | |
298b | 5.993 | 2550 | 4.17 | ||
Model4 | 300 | 5.945 | 2600 | 4.21 | |
Model30 | 300 | 5.944 | 2605 | 4.05 | |
Expt.2,a | 300 | 5.944 | 2610 | 4.06 | |
Expt.1 | 300 | 5.951 | 2608 | 4.46 | |
Expt.14 | 0 | 2653 | |||
Expt.31 | 2610 | 4.65c/3.95d |
These data are used to calculate the equations of state, which we report in Figs. S2 and S3 of the supplementary material.
Values estimated using a Debye model.
Ultrasonic experiment.
Shock wave experiment.
In Table II, we report the values of the ECs C11, C12, and C44 calculated with the three functionals together with the values of the bulk modulus, Young’s modulus, shear modulus, and Poisson’s ratio of polycrystalline molybdenum calculated using the Voigt–Reuss–Hill approximation. The temperature-dependent ECs have been measured in Refs. 14–17. Although there is not perfect agreement among these data, the 0 K values of Refs. 16–18 are quite close to each other. Taking as a reference the values of Ref. 18 extrapolated to 0 K, by adding the theoretical difference between 0 and 300 K, we find that the LDA errors for C11, C12, and C44 are 399 kbar (8%), 219 kbar (14%), and −27 kbar (−2%), while the PBE errors are −147 kbar (−3%), −7 kbar (−0.4%), and −105 kbar (−9%). PBEsol has errors 192 kbar (4%), 131 kbar (8%), and −40 kbar (−4%) smaller than LDA, but bigger than PBE. The PBE values of C11 and C12, and hence of its bulk modulus, are the closest to experiment. Since the calculation of the TDEC is computationally heavy, we calculated them using only this functional. Actually, as shown in Ref. 26, for several metals and as we confirmed in a recent study of tungsten,27 different functionals give different 0 K values of the ECs but the temperature and pressure dependence is almost independent of the functional.
The 0 K elastic constants calculated with the different functionals compared with experiment and one previous calculation. B, E, G, and ν are the bulk modulus, Young’s modulus, the shear modulus, and Poisson’s ratio, respectively.
. | T (K) . | a0 (a.u.) . | C11 (kbar) . | C12 (kbar) . | C44 (kbar) . | B (kbar) . | E (kbar) . | G (kbar) . | ν . |
---|---|---|---|---|---|---|---|---|---|
LDA | 0 | 5.884 | 5183 | 1815 | 1094 | 2938 | 3402 | 1301 | 0.307 |
PBEsol | 0 | 5.916 | 4976 | 1727 | 1081 | 2810 | 3318 | 1273 | 0.303 |
PBE | 0 | 5.974 | 4637 | 1589 | 1016 | 2605 | 3111 | 1196 | 0.301 |
PW919 | 0 | 5.988 | 4723 | 1604 | 1060 | 2644 | 3211 | 1237 | 0.297 |
Expt.14 | 0 | 4500.2 | 1729.2 | 1250.3 | 2653 | 3358 | 1303 | 0.289 | |
Expt.16 | 273.15 | 4637 | 1578 | 1092 | 2598 | 3232 | 1250 | 0.293 | |
Expt.16 (extrapolated) | 0 | 4800 | 1558 | 1124 | 2639 | 3354 | 1302 | 0.288 | |
Expt.15 | 300 | 4696 | 1676 | 1068 | 2683 | 3194 | 1227 | 0.302 | |
Expt.15 (extrapolated) | 0 | 4832 | 1656 | 1100 | 2715 | 3306 | 1275 | 0.297 | |
Expt.18 | 300 | 4648 | 1616 | 1089 | 2627 | 3222 | 1244 | 0.296 | |
Expt.18 (extrapolated) | 0 | 4784 | 1596 | 1121 | 2659 | 3334 | 1291 | 0.291 | |
Expt.19 | 300 | 2607 | 1251 |
. | T (K) . | a0 (a.u.) . | C11 (kbar) . | C12 (kbar) . | C44 (kbar) . | B (kbar) . | E (kbar) . | G (kbar) . | ν . |
---|---|---|---|---|---|---|---|---|---|
LDA | 0 | 5.884 | 5183 | 1815 | 1094 | 2938 | 3402 | 1301 | 0.307 |
PBEsol | 0 | 5.916 | 4976 | 1727 | 1081 | 2810 | 3318 | 1273 | 0.303 |
PBE | 0 | 5.974 | 4637 | 1589 | 1016 | 2605 | 3111 | 1196 | 0.301 |
PW919 | 0 | 5.988 | 4723 | 1604 | 1060 | 2644 | 3211 | 1237 | 0.297 |
Expt.14 | 0 | 4500.2 | 1729.2 | 1250.3 | 2653 | 3358 | 1303 | 0.289 | |
Expt.16 | 273.15 | 4637 | 1578 | 1092 | 2598 | 3232 | 1250 | 0.293 | |
Expt.16 (extrapolated) | 0 | 4800 | 1558 | 1124 | 2639 | 3354 | 1302 | 0.288 | |
Expt.15 | 300 | 4696 | 1676 | 1068 | 2683 | 3194 | 1227 | 0.302 | |
Expt.15 (extrapolated) | 0 | 4832 | 1656 | 1100 | 2715 | 3306 | 1275 | 0.297 | |
Expt.18 | 300 | 4648 | 1616 | 1089 | 2627 | 3222 | 1244 | 0.296 | |
Expt.18 (extrapolated) | 0 | 4784 | 1596 | 1121 | 2659 | 3334 | 1291 | 0.291 | |
Expt.19 | 300 | 2607 | 1251 |
In Fig. 1, we present the pressure-dependent ECs at 0 K calculated with the three functionals and compare them with the PBE results of Koči et al.21 There is a reasonable agreement between the two calculations, especially at low pressures. At 3000 kbar, our ECs are smaller than those of Koči et al. but quite close to them. At zero pressure, the pressure derivatives of the ECs are , , and , for all three functionals to be compared to the following experimental values:18 , , and .
Elastic constants as a function of pressure calculated within LDA (red lines), PBEsol (green lines), and PBE (blue lines) compared with the PBE results of Ref. 21.
Elastic constants as a function of pressure calculated within LDA (red lines), PBEsol (green lines), and PBE (blue lines) compared with the PBE results of Ref. 21.
We show in Fig. 2 the QHA isothermal and adiabatic ECs compared with the adiabatic experimental values. Keeping into account the zero point motion on both the lattice constant and the ECs themselves, we find C11 = 4564 kbar, C12 = 1589 kbar, and C44 = 996 kbar at 4 K, while computing the ECs from the strain derivatives of the energies at the lattice constant expanded by zero point motion effects within the QSA, we get C11 = 4593 kbar, C12 = 1565 kbar, and C44 = 996 kbar.
Quasi-harmonic isothermal (dashed lines) and adiabatic (solid line) elastic constants C11, C12, and C44 as a function of temperature compared with experimental adiabatic data from Ref. 14 (yellow circles), Ref. 16 (green circles), Ref. 15 (blue circles), and Ref. 17 (red circles).
As can be seen from Fig. 2, there is a good agreement between our calculated temperature dependence and the experimental data. From 24 to 2022 K, the experimental values17 decrease by 1264 kbar (27%), −63 kbar (−4%), and 231 kbar (21%) for C11, C12, and C44, respectively, while our values decrease by 1400 kbar (31%), −141 kbar (−9%), and 287 kbar (29%). In particular, the increase of C12 with temperature is found also in our QHA calculation, slightly overestimated with respect to experiment.
For comparison, we show in Fig. 3 the ECs calculated with the PBE functional within the QSA, which are in good agreement with those calculated in Ref. 20. In this case, from 4 to 2000 K, the decreases of C11, C12, and C44 are 380 kbar (8%), 127 kbar (8%), and 132 kbar (13%), respectively. The decrease of C11 and C44 is much smaller than in experiments (and within QHA), while C12 decreases with temperature instead of increasing as in experiment. We can understand this behavior using the QHA ECs calculated at a fixed volume, which do not contain any thermal expansion effect. For these ECs, C11 and C44 decrease with temperature, while C12 increases. Since QSA has only the effect of thermal expansion for C11 and C44, it misses the QHA contribution that gives a larger decrease, while for C12, it has no increasing QHA term. The QHA predicts an almost constant C12, which is the result of the cancellation between the decrease due to thermal expansion and the increase due to the use of the free energy derivatives instead of the energy derivatives.
Quasi-static isothermal (dashed lines) and adiabatic (solid line) elastic constants C11, C12, and C44 as a function of temperature compared with adiabatic experimental data from Ref. 14 (gold circles), Ref. 16 (green circles), Ref. 15 (blue circles), and Ref. 17 (red circles).
Using the QHA ECs, we have calculated the properties of polycrystalline molybdenum. In Fig. 4, we show the pressure dependence of the bulk modulus and the shear modulus in the range of pressures (up to 140 kbar) measured in Ref. 19. In addition to the 300 K calculation (green line), which can be compared with experiment, we show our predictions for 4, 1000, 1500, and 2000 K. We can see that the derivatives of the bulk and shear modulus with respect to pressure are well followed by our curves. Our values at 300 K are and against experimental values (obtained by a linear fit) and , respectively. These data are in agreement with the experimental values reported in Ref. 18: and and with the 0 K PBE theoretical results of Ref. 19: and . Regarding the temperature dependence of the adiabatic bulk and shear modulus, we find the following derivatives at 298 K: kbar/K and kbar/K (Fig. 4).
Adiabatic bulk and shear moduli of polycrystalline molybdenum against pressure computed at 5 K (red line), 300 K (green line), 1000 K (blue line), 1500 K (yellow line), and 2000 K (pink line), compared with room temperature experimental values of Ref. 19 (red circles).
Adiabatic bulk and shear moduli of polycrystalline molybdenum against pressure computed at 5 K (red line), 300 K (green line), 1000 K (blue line), 1500 K (yellow line), and 2000 K (pink line), compared with room temperature experimental values of Ref. 19 (red circles).
Finally, using Eqs. (9) and (10), we computed the compressional and shear sound velocities as a function of pressure for the same set of temperatures used in the previous picture. They are presented in Fig. 5. Even in this case, the pressure dependence of the sound velocity at 300 K is well reproduced by the calculation and the other curves are our prediction.
Compressional and shear sound velocities of polycrystalline molybdenum against pressure at 5 K (red line), 300 K (green line), 1000 K (blue line), 1500 K (yellow line), and 2000 K (pink line), compared with room temperature experimental values of Ref. 19 (red circles).
Compressional and shear sound velocities of polycrystalline molybdenum against pressure at 5 K (red line), 300 K (green line), 1000 K (blue line), 1500 K (yellow line), and 2000 K (pink line), compared with room temperature experimental values of Ref. 19 (red circles).
IV. CONCLUSIONS
We presented the temperature- and pressure-dependent thermoelastic properties of molybdenum calculated by the thermo_pw software and PAW pseudopotentials. We find that the QHA predicts the temperature dependence of the ECs in much better agreement with experiments than the QSA. Furthermore, we have used the QHA ECs to compute the pressure-dependent compressional and shear sound velocities in polycrystalline molybdenum (and the corresponding bulk and shear moduli). In addition to the calculation at 300 K, which is in good agreement with the experimental results of Liu et al.,19 we have calculated the low temperature (4 K) and high temperature (1000, 1500, and 2000 K) pressure-dependent curves, hoping that these calculations will stimulate an experimental investigation of these quantities.
For the sake of completeness, the phonon dispersions, the p–V equation of state at 300 and 2000 K, the temperature-dependent volume thermal expansion, the isobaric heat capacity, the adiabatic bulk modulus, and the average Grüneisen parameter for 0, 1000, 2000, and 3000 kbar have been calculated by the LDA, PBEsol, and PBE functionals, but since they are already available in the literature, we have moved them to the supplementary material.
SUPPLEMENTARY MATERIAL
The supplementary material that presents the calculated thermodynamic properties and a few tests on some numerical parameters is available at the following link.
ACKNOWLEDGMENTS
This work was supported by the Italian MUR (Ministry of University and Research) through the National Centre for HPC, Big Data, and Quantum Computing (Grant No. CN00000013). Computational facilities have been provided by SISSA through its Linux Cluster, ITCS, and the SISSA-CINECA 2021-2024 Agreement. Partial support has been received from the European Union through the MAX “MAterials design at the eXascale” Centre of Excellence for Supercomputing applications [Grant Agreement No. 101093374, co-funded by the European High Performance Computing Joint Undertaking (JU) and participating countries 824143]. X. Gong acknowledges the support received in the framework of the Joint Research Agreement for Magnetic Confinement Fusion between Eni and CNR.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
X. Gong: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). A. Dal Corso: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal).
DATA AVAILABILITY
The data that support the findings of this study are available within the article and its supplementary material.