Ab initio inspection of thermophysical experiments for molybdenum near melting

Refractory metals are of great importance because of their unique and desirable properties. First of all it concerns their high mechanical strength, high melting point, supreme resistance to heat and corrosion, and good electrical conductivity. Molybdenum is a refractory metal that plays an important role in metallurgy being employed to harden alloys used by aerospace and nuclear industries. It is also used as a pressure standard for ultra-high pressure static experiments.^1–3 Because of its essential role in the field of material and condensed matter science, molybdenum has attracted tremendous experimental and theoretical interest in its wide range of properties recently.

The main difficulty for experimental study of molybdenum in the liquid state is its high melting temperature (2.896 kK).⁴ Investigations of thermophysical properties of solid and liquid substances by conventional steady-state or quasi steady-state techniques are limited to below approximately 2500 K. Long exposure duration of a specimen at higher temperatures are limited due to increasing influence of associated processes such as chemical reactions, contamination, heat transfer, evaporation, loss of mechanical strength, electric insulation, etc.⁵ 

In this regard two main approaches that overcome some of these problems and thus extend the upper temperature limit for thermophysical measurements should be mentioned. The first one is a levitation technique, in which an external force (e.g. aerodynamic, acoustic, electromagnetic, electrostatic) is used to hold a material in space without a crucible.⁶ The absence of a crucible allows to avoid problems of the sample material interaction with its environment.⁷ Commonly used associated heating methods consist of inductive heating in case of an electromagnetic levitation, and laser heating for aerodynamic and electrostatic levitations. The levitation technique is very useful since it can handle liquid samples in containerless conditions and allow to measure a variety of thermophysical properties, such as heat capacity at constant pressure, density, surface tension, viscosity, etc.^8,9 However, it is limited to about 3 kK due to radiative heat losses. Moreover, the method can be sensitive to the chemical environment, and problems of stability and shape deformation of the levitated droplet can occur.⁵ 

The second technique is fast resistive pulse heating, in which a sample (wire, strips, etc.) is uniformly self-heated by electric current for a short period of time by the Joule effect. This approach minimises effects of thermal and chemical reactions with the surrounding medium and allows to reach much higher temperature ranges (up to 10 kK). Such dynamic heating technique can be classified according to various heating durations of the experiment and/or heating rates: millisecond, microsecond and submicrosecond (“wire explosion experiments”). The main advantages of using short duration experiments is preventing a metal rod from collapsing even after melting (temperature can rise up to 10,000 and even above for submicrosecond resolutions systems) and reducing of heat losses. However, the shorter pulse may lead to the lost of thermodynamic equilibrium and non-uniform heating due to the skin effect. An amount of thermophysical data that can be obtained using this method is impressive, and includes enthalpy of melting, specific heat at constant pressure, thermal expansion coefficient, sound speed, compressibility, transport properties (electrical and thermal conductivity, diffusivity, viscosity), critical data for the liquid–gas transition (in some cases), etc.¹⁰ 

Molybdenum has been studied extensively using the levitation technique and even more widely by the pulse heating method. However, its thermophysical properties in the vicinity of melting, especially in the liquid phase, are still not clear, because of strong contradictions in results of measurements. For example, density of molten molybdenum varies from 9.35 to 9.0 g/cm³ in dynamic experiments. The slope of the thermal expansion curve is also conflicting (from −0.44 to −0.75 g/(cm³ K)). Thus, the selection of experimental data is a main issue for constructing of a reliable equation of state of molybdenum in the lower-density and moderate pressure domain.¹¹ It also leads to the enormous scatter of estimates of the critical point — more than twofold on temperature (from 8 to 19.7 kK).

On the other hand, theoretical description of molybdenum in a condensed state is very challenging as well. Like many of the transition metals, it has a complex electronic structure that leads to the variety of unusual physical effects, like a positive sign of Seebeck and Hall coefficient near the ambient conditions, anomalous self-diffusion behaviour near the melting point,¹² dynamical stabilization of the crystal structure at high temperature by anharmonic effects,¹³ etc. This fact makes it virtually impossible to create an adequate analytic model of thermophysical properties of hot expanded molybdenum.¹⁴ 

Thus, it seems that recently the only available theoretical approach that can provide information about thermophysical properties of a substance in the region of hot expanded liquid and has a sufficient predictive power is ab initio or quantum molecular dynamics (QMD). This approach does not use any empirical data except for fundamental physical constants, so it can be used as a kind of a reference method for the analysis of available experimental data. Previously, this method was successfully implemented for the description of thermodynamical and transport properties of hot expanded nickel,¹⁵ aluminum,¹⁶ boron,¹⁷ calculation of the phase diagram of tantalum¹⁸ and estimation of tungsten critical point.¹⁹ 

In this paper we employ QMD simulations to obtain thermophysical properties of Mo in the region of hot expanded solid and liquid molybdenum. A detailed analysis of available experimental isobaric expansion data and extensive comparison over a wide set of properties with results of first-principle calculations is presented, possible reasons of contradictions are discussed.

We use method of quantum molecular dynamics (QMD) employing finite-temperature density functional theory (FT-DFT)²⁰ to obtain thermophysical properties of molybdenum. QMD simulations are performed in the Born-Oppenheimer approximation within the framework of the projector augmented-wave²¹ (PAW) method as implemented in the Vienna ab initio Simulation Package (VASP).^22–25 

The choice of the exchange-correlation (XC) functional is the most important aspect of DFT calculations,²⁶ that has the most significant impact on the result.²⁷ The generalized gradient approximation (GGA) of Perdew, Burke, and Ernzerhof^28,29 (PBE) is the most commonly used in condensed matter calculations, recently.³⁰ It should be also noted, that the PBE approximation provides good atomization energies, a property that is important for expanded states.³¹ For this reason we use GGA PBE for the XC functional in all simulations. However, it is known, that GGA systematically overestimates lattice constants. We calculated an equilibrium density at T = 300 K within the quasiharmonic approximation (QHA). The PHONOPY code³² was employed for performing QHA calculations. The equilibrium density at normal condition turned out to be $ρ0GGA=10.126$ g/cm³, that is 1% less than the experimental value⁴ of $ρ0exp=10.22$ g/cm³. The use of local density approximation (LDA) with the Ceperley-Alder (CA) parametrization³³ gave an even worse result: $ρ0LDA=10.58$ g/cm³, 3.5% more than the experimental value.

To reduce numerical efforts, the PAW potential with six valence electrons is used. The plane-wave cutoff energy is fixed at 400 eV.

All QMD simulations are made in the canonical (NVT) ensemble using the Nosé-Hoover³⁴ thermostat. The temperature of the electrons is fixed by the Fermi–Dirac distribution which defines the occupation numbers of bands. A sufficient number of bands is taken so that the occupation of the highest band was less than 10⁻⁵. Up to 800 electronic bands for calculations with 128 atoms are required. The time step is 2.0 fs.

The key aspect of accurate QMD calculation is numerical convergence, especially on the number of particles and k-points in the Brillouin zone. The convergence of the calculations of pressure and energy for the solid and liquid states at different densities and temperatures was extensively tested by varying the k-point sampling and number of atoms. We found out that a single k-point sampling is not sufficient for adequate description of thermophysical properties in a crystalline state. Based on our convergence tests we use the following settings: 128 atoms in the supercell and 2 × 2 × 2 Monkhorst-Pack³⁵ k-point set in the Brillouin zone for the solid state, and 128 atoms with the Baldereschi mean-value³⁶ k-point {1/4, 1/4, 1/4} for the liquid one. We collected representative examples of the convergence tests in Table I of  Appendix A.

Initial configuration for the calculations of the solid phase is the ideal bcc lattice, while the liquid phase simulations start from a disordered state. We do not determine melting temperature in this paper and use the reference value of T_m = 2.896 kK.⁴ This issue was studied by ab initio methods extensively,^37–40 and DFT is expected to be reliable.⁴¹ 

To estimate statistical uncertainty in the averaging of thermodynamic properties we performed a long QMD run for 128 atoms of Mo in the liquid state at a moderate rarefaction. As can be seen from Fig. 1(a), pressure fluctuates strongly during the simulation. Even the moving average over 1000 time steps varies within several kbars. Estimated standard deviations of the calculated average versus the trajectory length are shown in Fig. 1(b). At the averaging stage, each thermodynamic state was simulated during no less than 5000 steps in this work. Thus the standard error of a calculated average pressure can be estimated as about 1 kbar. The standard error of average energy estimated in a similar fashion is 0.3 kJ/mol.

The low-density domain of the phase diagram of molybdenum is extensively studied in experiments up to a temperature of 10 kK. The main source of experimental data in the liquid phase is the dynamic pulse heating method. A typical experiment consists of the rapid resistive self-heating of an electrical conductor specimen (typically a thin wire or rod of a ≈0.1–1 mm in diameter, and a few cm long) by the electrical current pulse. The base quantities measured in this type of experiments are voltage drop, current intensity, specimen diameter, surface radiation, and some others. The calculated quantities derived from the base ones consist of enthalpy vs. temperature, electrical resistivity vs. temperature, volume expansion vs. temperature or enthalpy, specific head capacity, thermal conductivity, and some others.⁵ An extensive ab initio study of thermoelectric transport properties of molybdenum is presented by French and Mattsson.²⁶ In this paper, we are going to focus on contradictions in density and enthalpy measurements.

To reduce the Rayleigh–Taylor instability in a liquid specimen the experimental chamber is usually filled with an inert gas at a static pressure of 0.15 GPa or greater.⁴² The measurement is essentially isobaric, corresponding to the initially established pressure.

Another type of experimental data which are available for liquid molybdenum in the vicinity of melting is non-contact thermophysical property measurements with an electrostatic levitation (ESL) technique.⁴³ This method uses the Coulomb force between a charged sample and electrodes and needs a high–speed feedback control system to stabilize the sample position. A high power Nd:YAG laser is used to melt the molybdenum sample. The metal properties are obtained in the conditions of high vacuum.

It was shown in Section II that absolute computational error might be commensurable with the experimental value of pressure due to the fluctuations of the pressure during the simulation. We can not significantly increase the number of atoms in the simulation to reduce oscillations, because the scaling of computational efforts is close to O(N³), where N is the number of valence electrons.^44,45 Therefore to reduce the influence of computational errors and to be able to reconstruct the isobaric expansion (IEX) curve at a given pressure with high accuracy we perform a series of calculations along isochors in the liquid phase and along the isotherm in the solid one. After that the calculated points are interpolated by a quadratic polynomial or linear fit, correspondingly, and interpolation relations are used for the 0 kbar isobar reconstruction or any other one if it is needed hereafter. Our calculations of isotherms (isochors) for the solid (liquid) state are presented in Fig. 2 (Fig. 3), correspondingly, and collected in Table II of  Appendix B.

Our calculations and available experimental data on isobaric expansion of molybdenum are presented in Fig. 4. Since the calculated density at normal conditions slightly differs from the experimental value, as it has been marked in Section II, it makes sense to compare the results of calculations and measurements in the units of relative density ρ/ρ₀, where ρ₀ is the density at normal conditions. In case of QMD $ρ0≡ρ0GGA=10.126$ g/cm³.

Solid state measurements are quite in agreement with each other. Pulse-heating data of Miiller and Cezairliyan⁴⁶ obtained by a millisecond resolution pulse-heating technique and Hixson and Winkler,⁴⁷ by sub-microsecond measurements, as well as ESL data of Paradis et al.^43,47 are presented. The QMD data are close to the Miiller and Cezairliyan approximation at 1.5 and 2.1 kK, however our calculations predict lower density of Mo at higher temperatures and agree almost perfectly with ESL measurements for the density of solid molybdenum in the range of 2.4–2.896 kK.

For the liquid state the experimental data from the following papers are shown: Seydel and Kitzel,⁴⁸ Hixson and Winkler,⁴⁷ Gathers,⁴² Pottlacher et al,⁴⁹ Hüpf et al,⁵⁰ Paradis et al.⁴³ These data are in a very strong contradiction as can be seen from Fig. 4. The density of liquid Mo observed in experiments is uniformly distributed from 9.35 to 9.0 g/cm³.

The best polynomial fit for our density data in the solid phase (1.5 ≤ T ≤ 2.896 kK) is

and in the liquid phase (2.896 ≤ T ≤ 6.2 kK) is

where T is in kK.

Surprisingly, QMD calculations agree excellently with the lowest curve of molten molybdenum density, presented in experimental works. Our points are close to the data by Pottlacher et al⁴⁹ and Hüpf et al.⁵⁰ The temperature dependence of liquid density is also in very good agreement with Refs. 49 and 50 up to 6 kK.

In experiment of Ref. 49 a static pressure of 2 kbar was applied and water was used as a pressure medium in a submicrosecond time resolution regime. The short circuit ringing period of the circuit was 6.2 μs, so heating rates greater than 10⁹ K/s were obtained. Some more detailed technical information is available in Ref. 51. In Ref. 50 experimental duration was about 50 μs corresponding to heating rates of 10⁸ K/s. Pressure was not specified.

Our calculations do not confirm anomalous behavior of the thermal expansion curve presented in Ref. 48. These data were obtained using a water-filled pressure vessel and short-circuited discharge time of 6 μs.

The question of the influence of water vapor propagating ahead of the expanding sample was raised in Refs. 42, 49, and 54. The reduced refractive index of the vapor could cause the apparent shadow size to increase. The application of a high pressure far above the critical point of water allows to reduce possible errors, as stated by Pottlacher et al.⁴⁹ 

The measurements by Hixson and Winkler⁴⁷ were made in a 2 kbar inert gas atmosphere. The diameter of the wire stock was 0.75 mm, that is more than in experiments of Refs. 48 and 49 (0.25 mm, correspondingly). While in the solid phase the agreement is good, the liquid–phase points of Hixson and Winkler are situated significantly higher.

The thermal expansion curves presented by Shaner et al.⁵² and Gathers⁴² have a steeper slope and predict higher densities of liquid molybdenum at the origin. Argon gas at a 2 kbar pressure was used in these experiments. The heating pulse length was tens of μs and the diameter of a wire was 1 mm.

It is interesting, that the data by Shaner et al.⁵² were revised later. Guillermet⁵³ took into account the systematic error in the measurements of temperature of liquid Ta reported by Gathers,⁵⁵ and used that information to estimate the corrections to the molar volume of liquid Mo assuming the error in the temperature measurements to be the same for both metals (the same apparatus was used). As can be seen from the figure, the corrected points are in much better agreement with our data and experiment by Pottlacher et al.⁴⁹ 

We do not reproduce the density of liquid Mo measured in ESL experiments⁴³ and obtained by Seydel and Kitzel.⁴⁸ One can distinguish two groups of data in the liquid phase in Fig. 4: with a moderate slope and with a big one. The first group contains our QMD points [slope −0.42 × 10⁻³ g/(cm³K)], experiments from Refs. 49 and 50 [−0.44 × 10⁻³ g/(cm³K)], and from Ref. 47 [−0.5 × 10⁻³ g/(cm³K)]. Experiments from Refs. 42 and 52 [−0.75 × 10⁻³ g/(cm³K)] and 48 [−0.78 × 10⁻³ g/(cm³K) can be assigned to the second group. The ESL curve⁴³ has an intermediate slope [−0.6 × 10⁻³ g/(cm³K)] between the two groups.

The detailed analysis of the sources of uncertainty of the density measurements in the levitation experiments was recently presented by Ozawa et al.⁵⁶ and Yoo et al.⁵⁷ The main difficulty concerns the estimation of an accurate droplet volume when the levitating droplet changes its form.⁵⁸ Translational oscillations of the droplet and high-temperature luminescence may also induce an incorrect estimation of the volume. For example density of liquid Zr measured with a levitation technique by separate authors differ by about 5%.⁵⁷ Refractory metals become luminous at elevated temperatures and it is hard to get an excellent contrast with the background. The effect of the background illumination on the measurement of the volume can be substantial.⁵⁹ The use of ultra-violet background illumination helps to reduce this effect as stated by authors⁶⁰ of Mo density measurements.

Complementary thermophysical data provided by ESL experiments are surface tension and viscosity, which are determined using the drop oscillation technique.^6,61 Actually, viscosity can be directly calculated by QMD by time integration of the auto-correlation function of the off-diagonal components of the stress tensor⁶² or from the transverse current correlation function.^63,64 It was also shown that the Stokes–Einstein (SE) relation for the Brownian motion is also valid for atoms:⁶⁵ 

and provides good agreement for the estimate of viscosity for liquid aluminum with direct methods and experiments.^62–64 Here k_B is the Boltzmann constant, D is a self-diffusion coefficient, and d is chosen to be the radius of the first peak in the radial distribution function. Coefficient 2 is used in Eq. (3) in the denominator instead of 6 as in the original SE relation.

Self-diffusion coefficient can be determined from the slope of the mean-square displacement at long times.

The use of the SE relation for our QMD data at T = 2.95 kK (d = 2.75 Å, D = 6.27 × 10⁻⁹ m²/s) gives a value of viscosity η ≈ 3.7 mPa s. It is lower than that measured by electrostatic levitation^66,67 (5 mPa s).

As was stated above, another thermophysical parameter available in dynamic heating experiments is enthalpy. Enthalpy difference under isobaric conditions of pulse heating experiments is equal to electrical power (the Joule effect). When temperature is measured as a function of time, one can relate the corresponding enthalpy measurements to obtain enthalpy as a function of temperature.^5,42

We can obtain the enthalpy change along the zero isobar ΔH ≡ H_T − H₃₀₀ = U_T − U₃₀₀, where U_T and U₃₀₀ are the internal energies of the system at a given temperature and at T = 300 K, correspondingly, directly from QMD calculations. To diminish the systematic error the convergence of U₃₀₀ was carefully tested (see Table I in  Appendix A).

Experimental and reference data^71–73 on specific enthalpy vs. temperature are presented in Fig. 5 together with the results of our QMD calculations.

The best polynomial fit for our enthalpy data in the solid phase (1.5 ≤ T ≤ 2.896 kK) is

and in the liquid phase (2.896 ≤ T ≤ 6.2 kK) is

where ΔH is in kJ/mol and T is in kK.

As can be seen from the figure, the discrepancy between the experimental data in the solid phase is negligible. QMD coincides excellently with the approximation from Standard Reference Materials report for SRM-781 molybdenum⁷⁰ up to T = 2.4 kK, however it provides slightly higher specific enthalpy values near the melting point. There is some contradiction between experimental data for enthalpy of liquid Mo. The calculated specific enthalpy of molten Mo at T_m is consistent with most experiments. Despite the excellent agreement in T—ρ plane (see Fig. 4), our data do not reproduce well the enthalpy measurements by Pottlacher et al.⁴⁹ However, quite good agreement is observed with data by Hixson and Winkler⁴⁷ at high temperatures. The slope of the linear regression fit of the QMD points in the range of 2.9–3.8 kK (or, in other words, at constant isobaric heat capacity C_p of liquid molybdenum 40.7 J/(mol K)) is in very good agreement with the reference value of IVTANTHERMO Handbook^71,72 (40 J/(mol K)) and recommended by Desai⁷⁴ (40.35 J/(mol K)). As can be seen from the figure, high temperature points significantly deviate from the linear regression. Thus, QMD predicts the increase of C_p of liquid Mo at high temperatures.

A rough estimate of C_p for solid Mo at the melting temperature (lower bound) made by linear regression of QMD data in the region 2.4–2.9 kK is 52 J/(mol K). It agrees well with IVTANTHERMO data (49.1 J/(mol K)) and the value recommended by Standard Reference Material report⁷⁰ (51.57 J/(mol K) at 2800 K). Recently, a comprehensive analysis of available experimental data on heat capacity of solid Mo was made by Bodryakov,⁷⁵ and the value 53.29 J/(mol K) was recommended.

The heat of fusion of molybdenum measured by pulse heating^{48,68,76–82} and by levitation calorimetry^69,83,84 methods varies in the range of 32.2—48.5 kJ/mol. We obtained a value of ΔH_fus = 34 kJ/mol and correspondingly entropy of fusion of ΔS_fus = 1.4 k_B/atom. Thus, our calculations confirm lower values of the enthalpy of melting for molybdenum.

Finally, data on specific enthalpy versus relative density is shown in Fig. 6. As can be expected, there is a strong contradiction between all experimental data in the liquid phase again, while there is good agreement in the solid one. We can not single out any experiment which our data reproduce in the liquid phase.

Since we determined enthalpy of fusion ΔH_fus and volume change ΔV_ls at melting it is natural to estimate the slope of the melting curve at P = 0 by the Clausius-Clapeyron relation:

We derive dT_m/dP = 47 K/GPa. Actually, the issue of the slope of molybdenum melting curve is a dramatic and highly debated^85–89 problem of a last few years. There is an enormous divergence in melting curves between laser-heated diamond anvil cell (DAC) experiments^90,91 and shock-wave (SW) measurements.^92–94 While DAC experiments showed nearly constant melting temperature as a function of pressure (not confirmed by theoretical studies^2,37,38,41), SW predicted a steep slope. In 2017, some light on the problem was shed by a novel DAC experimental work of Hrubiak et al.,⁹⁵ who presented a high-slope melting curve and suggested some possible explanations for the anomalously low melting observations. However, new results still point to a lower meting temperatures than those predicted by ab initio calculations.^37,38,41,96 Our result for dT_m/dP is higher than that predicted by Cazorla et al.⁴¹ (33.7 K/GPa) and Belonoshko et al. (2008)³⁸ (38 K/GPa). Better agreement is found with an earlier study by Belonoshko et al. (2004)³⁷ (40.5 K/GPa). It is much higher than the estimate based on DAC experiments by Errandonea⁹⁷ (10 K/GPa) and Hrubiak et al.⁹⁵ (21 K/GPa). However, it should be noted that in papers cited above high pressure melting points were fitted by the Simon functional form (a(1 + bP)^c) to get the slope at P = 0. The main and direct source of information about the slope at the origin of a melting curve is isobaric expansion experiments, mentioned and discussed earlier in this work. For example, an estimate of Shaner et al.⁵² is 33.3 K/GPa.

As it have been shown in this work, molybdenum is still one of the most complex metals for experimental and theoretical study, its thermophysical properties in vicinity of melting are highly ambiguous due to strong contradictions in available experimental data. Nevertheless, we have demonstrated that some experimental data can be accurately described by the QMD approach. In the solid phase the results of the first-principle calculations agree perfectly with the density measurements of Mo by the electrostatic levitation technique, while for liquid Mo excellent agreement is achieved with the measurements by the microsecond pulse heating technique by the group of the University of Technology of Graz (Pottlacher et al., Hüpf et al.). Previously we have demonstrated a similar situation for the density of liquid tungsten,¹⁹ the best agreement of QMD data was observed for the measurements of the workgroup of subseconds thermophysics in Graz.

Thus, our QMD calculations predict a substantial volume change on melting of about 5.5%, in addition to quite a low value of calculated enthalpy of fusion (34 kJ/mol). This fact leads to the steep slope of the melting curve at P = 0 predicted by QMD contrary to DAC experiments. We have obtained good agreement between the calculated isobaric heat-capacity of solid and liquid Mo near melting and experimental data. Moreover, we have observed the increase of C_p of liquid Mo at high temperatures.

Our findings are of significant importance for the construction of an accurate EOS for Mo in the vicinity of melting and in the subcritical region.

In conclusion we should say, that despite the difficulties in the interpretation of pulse heating experiments because of the complexity of concominant physical phenomena, problems of diagnostics and accurate pyrometry measurements,⁵ it is a unique technique that allows to obtain and measure a state of metals at kilokelvins.

We believe, that a combination of dynamic pulse heating and levitation experimental techniques, being accompanied by the QMD calculations will make it possible to obtain a complete set of accurate thermophysical data at the region of hot expanded substance and contribute significantly to the understanding of the structure and behaviour of matter under extreme conditions of temperature and pressure.

The majority of computations, development of codes, and treatment of results were carried out in the Joint Institute for High Temperatures RAS under financial support of the Russian Science Foundation (Grant No. 18-79-00346). The authors acknowledge the JIHT RAS Supercomputer Centre, the Joint Supercomputer Centre of the Russian Academy of Sciences, and the Shared Resource Centre “Far Eastern Computing Resource” IACP FEB RAS for providing computing time. Some numerical calculations were performed free of charge on supercomputers of Moscow Institute of Physics and Technology.

Table I contains representative examples of convergence tests. The standard error is estimated as 0.1 GPa for pressure, and 0.003 kJ/g for internal energy. The lattice constant obtained by QHA calculations at T = 300 K is used for QMD simulations at normal conditions. As can be seen, the convergence for energy is reached, while pressure is close to zero. We used the value U₃₀₀ = −10.933 kJ/g for ΔH calculation along the zero isobar. In the solid phase, the tests for isotherm T = 2.9 kK show that the simulations with 128 atoms and the Baldereschi’s mean-value k-point sampling³⁶ systematically underestimate the pressure and the energy values. More atoms in the supercell are necessary to reduce the systematic error. However, it turns out to be more computationally efficient to use 2 × 2 × 2 Monkhorst–Pack grid³⁵ for the Brillouin zone sampling as the number of symmetrically irreducible k-points is only 4, in this case. We checked that the application of the Γ-point-centered grid or a denser grid gives the same results. In the liquid phase, the tests for the isochors ρ = 8.875 g/cm³ and ρ = 8.75 g/cm³ demonstrate that the use of the Baldereschi’s mean-value point provides a good accuracy of thermodynamic properties for liquid molybdenum.

The EOS data derived from the QMD simulations are presented in Table II.