Ab initio simulation of liquid Mo and W near the liquid–gas coexistence curve

We present quantum molecular dynamics calculations of thermodynamic properties of expanded liquid molybdenum and tungsten. The liquid–gas coexistence curve and critical point parameters for these metals are estimated. We also analyze all available isobaric expansion experimental data and theoretical estimations of critical points for Mo and W.


Introduction
Thermodynamic properties of matter are described by an equation of state (EOS). The EOS is of considerable interest for basic research and has numerous important applications [1][2][3]. States of matter in the form of hot expanded liquid are characterized by high energy density and occupy a broad region in the phase diagram [4][5][6][7]. Our knowledge of these states is limited as both theoretical and experimental study of liquid metals at high temperature is hindered. Traditional methods of static thermophysical experiments cover the region of several thousands of Kelvins and several hundreds of bars. Investigation of thermodynamic properties of melted conductors in a wider region of temperatures, densities and pressures near the liquid-gas coexistence curve is possible by the isobaric expansion (IEX) technique. In this method a metal wire or stripe is rapidly heated by a powerful pulsed current and then expands into some medium, for example an inert gas or water. However the interpretation of IEX experiments is complicated because of the complexity of underlying physical phenomena [8,9]. This fact can explain a frequently occurring discrepancy between IEX experimental data for some metals, what is particularly vividly observed for refractory ones [2,10,11]. The determination of the critical point is also very important for the construction of the theory of liquid state. However direct measurements of critical parameters are possible only for low-melt metals (alkali or mercury) and only theoretical estimates can be done for the rest [12,13].
In this situation we hope that these issues can be clarified through the use of first-principles (ab initio) approaches. The most promising method for this task is quantum molecular dynamics (QMD) [14,15]. This approach can be considered as ab initio, because it does not employ any experimental data except for fundamental physical constants. Recent calculations of thermodynamic properties of liquid metals at atmospheric pressure demonstrated excellent agreement with experimental data [16]. So in this paper we are going to use QMD to describe thermal expansion of liquid molybdenum and tungsten up to the critical density. These metals are of particular interest, because some of available IEX experimental data contradict the others and existing estimates of the critical temperature differ more than twofold.

Methods and parameters of calculations
We use the QMD approach by employing the Vienna Ab initio Simulation Package (VASP) [17][18][19][20] to determine thermodynamic properties of liquid molybdenum and tungsten. For each QMD step the electronic density of valence electrons is calculated within the framework of finite-temperature density functional theory (DFT) [21], while the inner ones are replaced by an invariable core, and the ions are moved according to the forces acting on them which are determined via the Hellmann-Feynman theorem. The generalized gradient approximation (GGA) with the Perdew-Burke-Ernzerhof (PBE) [22,23] corrections for the exchangecorrelation functional was used. All calculations in this work were performed with PAW pseudopotentials [24]; 6 valence electrons are taken into account for both molybdenum and tungsten. The cutoff energy for the plane-wave basis was equal to 350 eV for both metals in majority of computations, but we increased this value to 450 eV at very low densities.
In most simulations of liquid molybdenum and tungsten 128 atoms were placed in a supercell (54 atoms at very low densities); the canonical ensemble with Nosé-Hoover thermostat was used. The evolution of the system was traced during 6 picoseconds. Thermodynamic quantities were taken as time averages of the molecular dynamics runs. The Brillouin zone was represented by the single Γ-point.
A method of successive approaching to the critical isotherm was used for the estimation of parameters of the critical point. The method is sketched in figure 1(a). We calculate a series of isotherms at gradually decreasing temperatures and analyze a density dependence of the derivative (∂P/∂ρ) T on each isotherm; here P -pressure, ρ-density, T -temperature. It is known that two conditions are valid on the critical isotherm at the critical point:   [25], open right triangles-Seydel and Kitzel [26], open circles-Pottlacher et al [27]; open triangles-Hixson and Winkler [28], dashed line-EOS by Fortov and Lomonosov [3]. Critical point estimates: red star-this work; filled square-Young and Alder [29]; filled pentagon-Fortov et al [30], filled diamond-Seydel and Fucke [31], filled down triangle-Seydel and Kitzel [26], filled triangle-Gathers [11], filled left triangle-Levashov et al [12], filled circle-Fortov and Lomonosov [3].
In a case of T < T c , the calculated point on an isotherm may occur inside a metastable or even thermodynamically unstable region. In this instance, the calculated isotherm has a local minimum and maximum. Since simulated QMD-systems are very small and evolve for only several picoseconds, instability effects do not have time to develop. Therefore we do not observe any clusters or cavities in the computational cell inside the thermodynamically unstable region.
Notice that the calculated critical pressure can be used for an estimation of the location of the liquid-gas coexistence curve, what is demonstrated in figure 1(b). Obviously, a subcritical isobar is closer to the binodal curve at high densities than a supercritical one, but the possibility of crossing the binodal should be carefully analyzed.

Results
We calculated a number of isotherms for expanded liquid molybdenum at densities from 7 to 1.25 g/cm 3 to determine the critical parameters. Our first estimations localized the search area in the range of 11-12 kK. The search for the critical isotherm of molybdenum turned out to be a quite challenging task because of the very narrow area of the characteristic change in the first pressure derivative on density. As can be seen from figure 2(a), QMD calculations predict the contact of the isotherm T c = 11.85 kK with the liquid-gas coexistence curve at ρ c = 1.6 g/cm 3 and P c = 7.1 kbar. This calculated critical pressure can be used to estimate the coexistence curve position.
The melting temperature of molybdenum is 2893 K at ambient pressure. This fact makes it difficult to obtain a reliable experimental information about properties of the liquid metal from static experiments. Molybdenum was studied using the IEX technique in which wire-shaped samples are resistively pulse heated in an inert gas atmosphere [25,27,28] or in water [26]. QMD isobars of molybdenum in comparison with the mentioned experimental data are presented in figure 2(b). As can be seen, the measurements by Shaner et al [25] and Pottlacher et al [27] [26]; open circles-Hixson and Winkler [34]; open left triangles-Koval et al [35]; dash-dotted line-soft-sphere approximation of data from [34] by Levashov et al [13]; dashed line-EOS by Fortov and Lomonosov [2]. Critical point estimates: red star-this work; filled square-Young and Alder [29]; filled pentagon-Fortov et al [30], filled down triangle-Seydel and Kitzel [26]; filled triangle-Fucke and Seydel [36]; filled left triangle-Levashov et al [13]; filled circle-Fortov and Lomonosov [2]. agree very well with each other, while the points by Seydel et al [26] demonstrate a steeper slope of thermal expansion in ρ-T diagram. Moreover an approximation curve of the data [26] has an anomalous positive second derivative. The data by Hixson et al [28] demonstrate higher densities of liquid molybdenum in the vicinity of the melting region in comparison with other experiments. Our QMD calculation of the zero-pressure isobar predicts the density of the molten molybdenum of about 9 g/cm 3 , that agree well with experiments of fast resistive pulse heating in inert gases. The slope of the calculated curve does not confirm the anomalous behavior of the thermal expansion curve predicted by Seydel et al [26] and demonstrate excellent agreement with experiments [25,27]. We also can conclude that our estimation of the coexistence curve of molybdenum differ from the existing Lomonosov's EOS [3].
The available critical point estimates for molybdenum spread from 8 to 16.5 kK [3,11,[29][30][31][32][33], however the discrepancy in the critical density is considerably smaller and varies within 2.3-3.2 g/cm 3 . The Lomonosov's EOS predicts a slightly higher value of ρ c = 3.7 g/cm 3 . Thus our estimation of the critical density of 1.6 g/cm 3 is lower than in other works. As can be seen from figure 2(b) the evaluation of the critical point from [31] demonstrates better agreement with results of our QMD calculations than the others.
The method of successive approaching to the critical isotherm by QMD calculations for liquid tungsten gives an estimation for the critical temperature as T c = 14 kK. The analysis of the first density derivative of this isotherm gives the parameters of the critical points ρ c = 2.9 g/cm 3 and P c = 10.4 kbar, see figure 3(a). The critical and zero pressure isobars of tungsten obtained by QMD calculations are presented in figure 3(b). Similar to molybdenum, it is not surprising that properties of liquid tungsten are studied only in dynamic experiments, because tungsten has the highest melting temperature among all metals. Available IEX data [26,34,35] are also presented in figure 3(b). As can be seen from the figure there is a significant contradiction between the IEX experimental data for liquid tungsten. The most considerable difference is observed between the 5 1234567890 ''""  [35] and earlier experiments of Seydel et al [26] and Hixson et al [34]. The density of molten tungsten given in [35] is almost 0.5 g/cm 3 higher than that in the previous experiments. The results of our QMD calculations of the zero pressure isobar are not consistent with the data from [35], but they are in excellent agreement with density of molten tungsten observed in the earlier experiments. The slope of the calculated thermal expansion curve is slightly softer than in [34]. Again, similar to molybdenum, there is no agreement between available estimations of the critical point for tungsten. Existing estimates varies from 13.4 to 21 kK [2,29,30,32,36]. Our estimation of the critical point parameters agrees with the prediction of Fucke and Seydel [36], however the value of the critical density is about 1.4 g/cm 3 lower. Nevertheless such a low value of the critical density agree well with the soft-sphere approximation of experimental data from [34] given by Levashov et al [13].

Conclusions
In this work, we have performed QMD calculations of thermodynamic properties of expanded liquid molybdenum and tungsten. The QMD-calculated zero pressure and critical isobars have been presented in the range of densities from 9 to 1.6 g/cm 3 for Mo and from 16.4 to 2.9 g/cm 3 for W. Our estimation of the critical point for Mo is T c = 11.85 kK, ρ c = 1.6 g/cm 3 and P c = 7.1 kbar, and for W is T c = 14 kK, ρ c = 2.9 g/cm 3 and P c = 10.4 kbar. We have made a comparison with all available isobaric expansion experimental data and theoretical estimations of the critical points for the metals under consideration. Our QMD calculations do not confirm the anomalous behavior of the thermal expansion curve of molybdenum predicted by Seydel and Kitzel [26]. The density of molten tungsten given by Koval et al [35] is also not consistent with the results of our simulations.