# Melting curve and Hugoniot of molybdenum up to 400 GPa by ab initio simulations

###### Abstract

We report ab initio calculations of the melting curve and Hugoniot of molybdenum for the pressure range GPa, using density functional theory (DFT) in the projector augmented wave (PAW) implementation. We use the “reference coexistence” technique to overcome uncertainties inherent in earlier DFT calculations of the melting curve of Mo. Our calculated melting curve agrees well with experiment at ambient pressure and is consistent with shock data at high pressure, but does not agree with the high pressure melting curve from static compression experiments. Our calculated and Hugoniot relations agree well with shock measurements. We use calculations of phonon dispersion relations as a function of pressure to eliminate some possible interpretations of the solid-solid phase transition observed in shock experiments on Mo.

## 1 Introduction

The melting curves of transition metals at pressures up to the megabar region are highly controversial, particularly for b.c.c. metals. Diamond anvil cell (DAC) measurements find that the melting temperature increases by only a few hundred K over the range 1 - 100 GPa [1, 2], while shock experiments indicate an increase of several thousand K over this range [3, 4, 5]. There have been several ab initio calculations on transition-metal (P) curves, and the predictions agree more closely with the shock data than with the DAC data [6, 7, 8, 9]. A challenging case is molybdenum, where there are very large differences between DAC and shock data [12], and where the shock data reveal two transitions, the one at high pressure ( 380 GPa) being attributed to melting, and the one at low pressure ( 210 GPa) to a transition from b.c.c. to an unidentified structure [4]. We report here on new ab initio calculations of () for Mo, and on the and relations on the Hugoniot. We also report preliminary information that may help in searching for the unidentified high-P solid phase of Mo.

We use density functional theory (DFT), which gives very accurate predictions for many properties of transition metals, including Hugoniot curves [10]. DFT molecular dynamics (m.d.) was first used to study solid-liquid equilibrium 12 years ago [11], and several different techniques are now available for using it to calculate melting curves. In such calculations, no empirical model is used to describe the interactions between the atoms, but instead the full electronic structure, and hence the total energy and the forces on the atoms, is recalculated at each time step. There have been earlier DFT calculations on the melting of Mo, but the techniques used were prone to superheating effects [6, 9]. In the present work, we use the so-called “reference coexistence ” techniques [19, 17, 14], which does not suffer from this problem.

Our work has several aims. First, we want to improve on the reliability and accuracy of the predicted melting curve of Mo obtained from DFT; second, we use DFT to predict the and relations on the shock Hugoniot; third, we want to identify the unknown solid phase of Mo observed in shock experiments. Our tests on the accuracy of DFT for Mo, and our extensive calculations of the Mo melting curve will be reported in detail elsewhere [18], so we present only a summary here. However, our very recent calculations on the shock Hugoniot will be presented in more detail. These are important, because temperature is very difficult to measure in shock experiments [35] and DFT gives a way of supplying what is missing in the shock data. Our search for the unidentified solid structure of Mo is at the exploratory stage, but we present results for phonon frequencies as a function of pressure, which allow us to rule out some possibilities.

In the following, we summarise our tests on the accuracy of the DFT techniques (Sec. 2), and outline the reference coexistence technique. In Sec. 3 we present our results for the DFT melting curve and Hugoniot of Mo up to 400 GPa, and the study of the phonon frequencies. Discussion and conclusions are in Sec. 4.

## 2 Techniques and tests

Our DFT calculations are performed mainly with the projector augmented wave (PAW) implementation [20], using the VASP code [21, 25], since PAW is known to be very accurate. The main uncontrollable approximation in DFT is the form adopted for the exchange-correlation functional . To test the accuracy of PAW, and the effect of , we have compared our predictions for the relation of the b.c.c. Mo crystal against low- experimental results (Fig. 1). The pressure predictions from PAW using the Perdew-Burke-Ernzerhof (PBE) [28] and local-density approximation (LDA) [27] forms of deviate by % in opposite directions from the experimental data, but we adopt the PBE form, because the deviations in this case are rather constant. The PBE results of Fig. 1 were obtained with 4s states and below in the core and all other states in the valence set. Inclusion of 4s states in the valence set makes no appreciable difference to the PBE results. We also tested the PAW implementation itself by repeating the calculations with the even more accurate full-potential linearized augmented plane-wave (FP-LAPW) technique [22, 23, 24], using the WIEN2k code [33]. As shown in Fig. 1, PAW and FP-LAPW results are almost identical. Further confirmation for the accuracy of the PAW implementation and the PBE functional comes from our comparisons of the calculated phonon dispersion relations for the ambient-pressure b.c.c. crystal with experiment (Fig. 2).

The “reference coexistence” technique for calculating ab initio melting curves has been described elsewhere [18], but we recall the main steps. First, an empirical reference model is fitted to DFT m.d. simulations of the solid and liquid at thermodynamic states close to the expected melting curve. Then, the reference model is used to perform simulations on large systems in which solid and liquid coexist, to obtain points on the melting curve of the model. In crucial third stage, differences between the reference and DFT total energy functions are used to correct the melting properties of the reference model, to obtain the ab initio melting curve. In the present work, the total energy function of the reference model is represented by the embedded-atom model (EAM) [30, 31], consisting of a repulsive inverse-power pair potential, and an embedding term describing the d-band bonding. The detailed procedure for fitting to DFT m.d. data will be reported elsewhere [18], but we note that we needed to re-fit the model in different pressure ranges.

The simulations on solid and liquid Mo in stable thermodynamic coexistence using the fitted reference model employed systems of 6750 atoms, and we checked that the results do not change if even larger systems are used. The protocols for preparing these simulated systems, and for achieving and monitoring stable coexistence were similar to those used in earlier work on the melting curve of Cu [17]. The procedures for correcting the melting curve of the reference model, which depend on calculations of the the free energy differences between the reference and DFT systems, have been described and validated in earlier work [18].

## 3 Results

### 3.1 Melting properties

Fig. 3 shows the melting curve of our reference EAM model, and the melting curve obtained from this by correcting for the differences between the DFT and reference total-energy functions; earlier ab initio-based calculations of the Mo melting curve due to Moriarty and to Belonoshko et al. are also indicated [6, 9]. We also show points on the melting curve from DAC and shock measurements [1, 4]. The differences between reference and corrected melting curves are only a few hundred K, so that the corrected curve should be very close to the melting curve that would be obtained from the (PBE) exchange-correlation functional if no statistical-mechanical approximations were made. Because we avoid approximations of earlier ab initio-based work, our present results should be a more accurate representation of the DFT melting curve. Up to 100 GPa, the differences between our results and earlier DFT work are rather small, and we confirm that DFT gives a much higher melting slope than that given by DAC experiments. We obtain an accurate value of at by fitting our melting curve to the Simon equation , with K, GPa, . The resulting value of K is close to the accepted experimental value K. Our value of K GPa at agrees with an older experimental value of K GPa [32]. Our DFT melting curve is consistent with the point obtained at GPa from shock measurements. However, we stress that the temperature of the “experimental” point was not measured, but estimated by considering superheating corrections to the shock-wave data [12, 13]. Because of this, it is important to try and corroborate the estimated experimental temperature, and this can be done by ab initio calculations, as we explain next.

### 3.2 Hugoniot curves

The pressure , volume and internal energy in the shocked state are related to the initial volume and internal energy by the well-known Rankine-Hugoniot formula [34]:

(1) |

Since the internal energy and pressure are given in terms of the Helmholtz free energy by and , we can calculate the Hugoniot from our DFT simulations, provided we can calculate as a function of and . So far, we have done this only for the b.c.c. crystal in the harmonic approximation, in which . Here, is the free energy of the rigid perfect crystal, including thermal electronic excitations, and is calculated from the phonon frequencies ( the wavevector, the branch). We calculate in the classical limit, in which per atom, with and is the geometric average of phonon frequencies over the Brillouin zone. The methods used for to calculate and the frequencies were similar to those used in our earlier work on Fe (Ref. [16]). For a set of temperatures, we calculated at a set of volumes, and fitted the volume dependence with a third-order Birch-Murnaghan equation [36]. The temperature dependence of the coefficients in this equation were then fitted with a third-order polynomial. The phonon frequencies were calculated at 12 volumes in the range Å/atom, as explained elsewhere [18]. The volume dependence of the average was then fitted with a third-order polynomial.

To obtain and from our fitted free energy, for each value of we seek the at which the Rankine-Hugoniot equation is satisfied, and from this we obtain . For and , we used values from our GGA calculations; we checked that use of the experimental made no significant difference. Comparison of our calculated with measurements of Hixson et al. [4, 5] (left panel of Fig. 4) shows excellent agreement. In the right panel, we compare our with both uncorrected results of Hixson et al. and also with results corrected for superheating, and our results confirm their temperature estimates.

### 3.3 Solid-solid phase transition

Efforts have been made to identify the high-/high- structure of Mo indicated by shock experiments [4]. Hixson et al. used their theoretical prediction of b.c.c. h.c.p. phase transition at GPa and K to suggest the h.c.p structure. However, later calculations locate this transition at higher pressures ( GPa), and temperature stabilisation of h.c.p. over b.c.c. below melting seems improbable. Furthermore, there are recent claims that under pressure b.c.c Mo transforms first to f.c.c. rather than h.c.p. [37]. Very recently, Belonoshko et al. [9] reported ab initio simulations on the f.c.c and A15 structures at temperatures and pressures near those of interest (namely, GPa and K), and concluded that both structures are unlikely high- phases of Mo. We have performed ab initio calculations similar to those of Belonoshko et al. on the h.c.p. and structures ( GPa and K). We find that temperature neither favours any of them respect to b.c.c., thus they must be excluded as well. Recently, a 2nd-order phase transition from cubic to rhombohedral has been observed in vanadium at GPa [38]. This appears to be related to an earlier ab initio prediction of a phonon softening in V [39]. This suggests that a similar structural transition might occur in Mo. We have used our calculated phonon dispersion relations of Mo over the range GPa to test this. Fig. 5 depicts our results at 0, 77, 136, 274 and 400 GPa, but we see that no phonon anomaly like that reported for V is observed at any wavevector. This indicates that we should rule out structures based on small distortion of b.c.c.

## 4 Discussion

An important outcome of the present work is improved DFT calculations of the melting curve of Mo over the pressure range GPa. In particular, our techniques allow us to avoid the superheating errors which appear to affect an independent recent DFT study of Mo melting [9]. The accuracy of our calculations is confirmed by the very close agreement with experiment for the melting temperature and the melting slope at ambient pressure. The results fully confirm that the increase of by K over the range GPa predicted by DFT is about 10 greater than that deduced from DAC measurements. A second important outcome is that our calculations of the temperature along the shock Hugoniot support earlier temperature estimates based on experimental data but corrected for possible superheating effects [12]. This allows us to compare more confidently the point on the melting curve at GPa derived from shock data with our predicted melting curve, and we confirm that at this pressure is ca. 8650 K. This is far above any reasonable extrapolation of the DAC data. Concerning the search for the unknown crystal structure of Mo indicated by shock experiments to exist above ca. 220 GPa, we have been able so far only to rule out some possibilities. Our calculations of the phonon dispersion realations in the b.c.c. structure over the range GPa reveal no softening of any phonons, and no indication of any elastic instability. This means that the new crystal structure does not arise from small distortion of b.c.c. This is interesting in the light of the recent discovery of the elastic instability of b.c.c. V above GPa, predicted initially by DFT, and observed very recently in x-ray diffraction experiments [38]. It seems that the structural transition in Mo is of a different kind.

The large conflict between the melting curve of Mo derived from DAC measurements on one side and from shock experiments and DFT calculations on the other side must be due either to a misinterpretation of the DAC data or to a combination of serious DFT errors and misinterpretation of shock data. Given the accuracy of DFT that we have been able to demonstrate (low-temperature curve, Hugoniot curve, ambient- phonons), we think there is little evidence for significant errors in DFT, which is also in good accord with shock data. A possible explanation might be that the large temperature gradients and non-hydrostatic stress in DAC experiments might give rise to flow of material giving the appearance of melting, even well below the thermodynamic melting temperature. We also note recent evidence that temperature measurement in DAC experiments may be subject to previously unsuspected errors [40], though probably not of the size needed to resolve the conflict by themselves.

The work was supported by EPSRC grant EP/C534360, which is 50% funded by DSTL(MOD), and by NERC grant NE/C51889X/1. The work was conducted as part of a EURYI scheme award to DA as provided by EPSRC (see www.esf.org/euryi).

## References

## References

- [1] Errandonea D., Schwager B., Ditz R., Gessmann C., Boehler R. and Ross M. 2001 Phys. Rev. B 63 132104
- [2] Errandonea D., Somayazulu M., Häusermann D. and Mao D. 2003 J. Phys: Condens. Matter 15 7635
- [3] Brown J. M. and Shaner J. W. Shock Waves in Condensed Matter 1983, ed. J. R. Assay, R. A. Graham and G. K. Straub (Amsterdam, Elsevier, 1983) p. 91
- [4] Hixson R. S., Boness D. A., Shaner J. W. and Moriarty J. A. 1989 Phys. Rev. Lett. 62 637
- [5] Hixson R. S. and Fritz F. N. 1992 J. Appl. Phys. 71 1721
- [6] Moriarty J. A. 1994 Phys. Rev. B 49 12431
- [7] Moriarty J. A., Belak J. F., Rudd R. E., Söderlind P., Streitz F. H. and Yang L. H. 2002 J. Phys: Condens. Matter 14 2825
- [8] Wang Y., Ahuja R. and Johansson B. 2002 Phys. Rev. B 65 014104
- [9] Belonoshko A. B., Simak S. I., Kochetov A. E., Johansson B., Burakovsky L. and Preston D. L. 2004 Phys. Rev. Lett. 92 195701
- [10] Gillan M. J., Alfè D., Brodholt J., Vočadlo and G. D. Price 2006 Rep. Prog. Phys. 69 2365
- [11] Morris J. R., Wang C. Z., Ho K. M. and Chan C. T. 1994 Phys. Rev. B 49 3109
- [12] Errandonea D. 2005 Physica B 357 356
- [13] Luo S. -N. and T. J. Ahrens 2003 Appl. Phys. Lett. 82 1836
- [14] Alfè D., Gillan M. J. and Price G. D. 1999 Nature 401 462
- [15] Vočadlo L. and Alfè D. 2002 Phys. Rev. B 65 214105
- [16] Alfè D., Price G. D. and Gillan M. J. 2001 Phys. Rev. B 64 045123
- [17] Vočadlo L., Alfè D., Price G. D. and Gillan M. J. 2004 J. Chem. Phys. 120, 2072
- [18] Cazorla C., Gillan M. J., Taioli S. and Alfè D. 2007 J. Chem. Phys., in press.
- [19] Alfè D., Gillan M. J. and Price G. D. 2002 J. Chem. Phys. 116 6170
- [20] Blöchl P. E. 1994 Phys. Rev. B 50 17953
- [21] Kresse G. and Joubert D. 1999 Phys. Rev. B 59 1758
- [22] Andersen O. K. 1975 Phys. Rev. B 12 3060
- [23] Koeling D. D. and Arbman G. O. 1975 Journal of Phys. F 5 2041
- [24] Singh D. 1994 Kluwer Academic Publishing ISBN 0-7923-9421-7
- [25] Kresse G. and Furthmüller J. 1996 Phys. Rev. B 54 11169
- [26] Zarestky J., Stassis C., Harmon B. N., Ho K. M. and Fu C. L. 1983 Phys. Rev. B 28 697
- [27] Ceperley D. M and Alder B. I. 1980 Phys. Rev. Lett. 45 566
- [28] Perdew J. P., Burke K. and Ernzerhof M. 1996 Phys. Rev. Lett. 77 3865
- [29] Verma A. K., Rao R. S. and Godwal B. K. 2004 J. Phys.:Condens. Matter 16 4799
- [30] Daw M. S. and Baskes M. I. 1984 Phys. Rev. B 29 6443
- [31] Finnis M. W. and Sinclair J. E. 1984 Phil. Mag. A 50 45
- [32] Shaner J. W., Gathers G. R. and Minichino C. 1977 High Temp. High Pressures 9 331
- [33] Blaha P., Schwarz K., Madsen G. K., Kvasnicka D. and Luitz J. 2001 WIEN2k: An Augmented Plane Wave plus Local Orbital Program for Calculating Crystal Properties Technical University of Vienna
- [34] Poirier J. P. in Introduction to the Physics of the Earth’s Interior (Cambridge University Press, Cambridge, England, 1991), see e.g., Chap. 4.
- [35] Yoo C. S., Holmes N. C., Ross M., Webb D. J. and Pike C. 1993 Phys. Rev. Lett. 70 3931
- [36] Birch F. 1978 J. Geophys. Res. 83 1257
- [37] Jona F. and Marcus P. M. 2005 J. Phys.: Condens. Matter 17 1049
- [38] Ding Y., Ahuja R., Shu J., Chow P., Luo W. and Mao H. 2007 Phys. Rev. Lett. 98 085502
- [39] Suzuki N. and Otani M. 2002 J. Phys.: Condens. Matter 14 10869
- [40] Benedetti L. R., Guignot N. and Farber D. L. 2007 J. of Appl. Phys. 101 013109