Crystal Structures and Exotic Behavior of Magnesium under Pressure

J. Phys. Chem. C 2010, 114, 21745–21749

21745

Crystal Structures and Exotic Behavior of Magnesium under Pressure Peifang Li,†,‡ Guoying Gao,† Yanchao Wang,† and Yanming Ma*,† State Key Lab of Superhard Materials, Jilin UniVersity, Changchun 130012, People’s Republic of China, and College of Physics and Electronic Information, Inner Mongolia UniVersity for the Nationalities, Tongliao 028043, People’s Republic of China ReceiVed: August 27, 2010; ReVised Manuscript ReceiVed: October 25, 2010

The high-pressure structures of magnesium (Mg) have been extensively explored through our newly developed particle swarm optimization algorithm on crystal structural prediction. Two structures with face-centered cubic (fcc) and simple hexagonal (sh) symmetries are discovered to be stable at 456-756 GPa and above 756 GPa, respectively. Especially, the sh structure, which is known to occur at high pressure in P, Si, and Ge, is reported here for Mg. This structure can be derived from the fcc lattice by distortion of the R and γ angles from 60° to 90°. More intriguingly, the calculated valence electron localization function reveals an electride nature of the fcc and sh structures with valence electrons localized in the interstitial regions, analogous to what was recently reported in Li, Na, K, and Ca under high pressure. However, what makes Mg unique is that it remains metallic. The temperature-pressure phase diagram of Mg has also been explored using quasiharmonic approximation, and the finite-temperature phase boundaries of the fcc and sh structures are determined. Introduction The high-pressure behavior of alkali and alkaline-earth metal elements has long attracted great attention. At ambient conditions, these metallic elements adopt simple structures, such as face-centered cubic (fcc), body-centered cubic (bcc), or hexagonal close-packed (hcp).1 Under pressure, they show extremely complex or even incommensurate structures2,3 with unusually large unit cells and low symmetry. At even higher pressures, the reappearance of close-packed structures was reported in Ba,4 Na,5 K,6 Rb,7 Cs,7 and Ca.8 With the change in crystal structure under pressure, alkali and alkaline-earth metal elements exhibit unusual physical properties, for example, strong departure from the free-electron-like behavior, valence bandwidth narrowing, decrease of coordination numbers, and increase of valence electron density in the interstitial regions.5,9 Mg adopts the hcp structure10 at ambient conditions. It transforms into bcc at 50 GPa and room temperature,11 but into a double hcp structure above 9.6 GPa at 1277 K.12 Above 70 GPa, Moriarty and McMaHan13 predicted that bcc might transform into fcc either at 180 GPa by generalized pseudopotential theory (GPT) or at 790 GPa by the linear-muffin-tin orbital (LMTO) method. The large discrepancy in the transition pressures between two different methods is quite surprising and remains a mystery. Later, Ahuja et al.14 predicted that the simple cubic (sc) structure might exist above 660 GPa. It is clear that the high-pressure structures of Mg are still elusive. It is also worth mentioning that the high-pressure behavior of Mg plays an important role in planetary physics, being that MgSiO3 is the planet-forming silicate stable at pressures and temperatures beyond those of Earth’s core-mantle boundary and dissociates into MgO and SiO2 at pressures (∼1000 GPa) and temperatures expected to occur in the cores of gas.15 The exploration of the crystal structures and physical properties of Mg under ultrahigh pressure is thus greatly desirable. * To whom correspondence should be addressed. E-mail: [email protected]. † Jilin University. ‡ Inner Mongolia University for the Nationalities.

In this paper, we employ our newly developed particle swarm optimization (PSO) technique on crystal structure prediction16 through the CALYPSO (Crystal structure AnaLYsis by Particle Swarm Optimization) code to extensively explore the crystal structures of Mg under ultrahigh pressure. We have revealed that fcc is the high-pressure structure of Mg and clarified that the bcc f fcc transition pressure is 456 GPa. At higher pressure, an sh phase is stable above 756 GPa. The two structures of Mg have deviated from “ideal” metallic behavior, and their valence electrons are mostly repulsion by core electrons into the lattice interstitial sites. We have constructed the phase diagram of Mg to 850 GPa, which could be greatly helpful for future experiments. Computational Approach PSO is a computational method that optimizes a problem by iteratively improving a candidate solution with regard to a given measure of quality. It is originally attributed to Kennedy, Eberhart, and Shi and was intended for simulating social behavior in the mid 1990s.17,18 A basic variant of the PSO algorithm works by using a population (called a swarm) of candidate solutions (called particles). These particles are selfadjusted in the search-space by the following operations t+1 t t+1 xi,j ) xi,j + Vi,j

(1)

t+1 t t t t t Vi,j ) ωVi,j + c1r1(pbesti,j - xi,j ) + c2r2(gbesti,j - xi,j ) (2)

The movements of the particles are guided by their own best known position in the search-space as well as the entire swarm’s best known position. PSO has been verified to perform well on many optimization problems.19,20 For the first time, we have implemented the PSO algorithm on crystal structure prediction in the CALYPSO code.16 Our global minimization method through the CALYPSO code for predicting crystal structures comprises mainly five features: (1) CALYPSO is based on an

10.1021/jp108136r  2010 American Chemical Society Published on Web 11/18/2010

21746

J. Phys. Chem. C, Vol. 114, No. 49, 2010

Li et al.

efficient global minimization of free energy surfaces merging total-energy calculations via the PSO technique; (2) CALYPSO requires only chemical compositions for a given compound to predict stable or metastable structures at given external conditions (e.g., pressure); (3) a particularly devised geometrical structure parameter was implemented to eliminate similar structures during structure evolution for enhancing the structure search efficiency; (4) a variable unit cell size technique was designed to greatly reduce the computational cost; and (5) CALYPSO imposes the symmetry constraint in the structure generation to realize diverse structures and reduce search space and optimization variables; thus, the global structure convergence can be significantly fastened. This methodology has been successfully applied to various known experimental structures on elemental, binary, and ternary compounds with metallic, ionic, and covalent bonding. It is proved to be powerful with high efficiency and a high success rate.16 Variable-cell PSO simulations with up to 16 atoms in the unit cell are performed for Mg at 0, 20, 200, 500, 800, and 1000 GPa. The underlying ab initio structural relaxations were performed using density functional theory within the PerdewBurke-Ernzerhof generalized gradient approximation,21 as implemented in VASP (Vinna Ab Initio Simulation Package) code.22 The all-electron projector-augmented wave (PAW) method23 was adopted with the PAW potentials taken from the VASP library24 where 2p63s2 states were used as valence electrons. The cutoff energy (600 eV) for the expansion of the wave function into plane waves and proper Monkhorst-Pack k meshes were chosen to ensure that all structures are well converged to better than 1 meV/atom. The results were also double-checked by full-potential linearized augmented planewave (FP-LAPW) methods within the WIEN2K codes.25 The temperature-pressure phase diagram has been explored using quasi-harmonic approximation. The Helmholtz energy F at volume V and temperature T can be approximated as

F(V, T) ) E(V) + Fvib(V, T)

(3)

Under quasi-harmonic approximation, Fvib can be calculated from phonon DOS by ∞

Fvib(V, T) ) NrkBT

[ ( )]

∫ dωg(ω, V)ln 2 sinh 2kpωBT 0

(4)

where F is the free energy of the primitive unit cell, N is the number of primitive unit cells, r is the number of degrees of freedom in the primitive unit cell, ω is the phonon frequency, g(ω) is the phonon DOS at frequency ω, p is the Planck constant, and kB is the Boltzmann’s constant. The phonon frequencies were calculated by the direct supercell approach, which uses the forces obtained by the Hellmann-Feynman theorem calculated from the optimized supercell. Results and Discussion We predicted correctly the hcp and bcc structures at 0 and 200 GPa, respectively. At 20 GPa, we also uncovered the double hcp structure proposed by Errandonea et al.;12 however, it is much higher in enthalpy than the hcp and bcc structures. At 500 GPa, our simulation predicted that the most stable structure is the fcc structure, as earlier proposed by Moriarty and McMaHan.13 At 800 GPa, an sh phase (space group P6/mmm,

Figure 1. (a) Enthalpy curves (relative to the bcc structure) as a function of pressure. Inset: enthalpy for the fcc structure relative to the bcc structure using the FP-LAPW method. (b) Theoretical equation of states of the hcp (dashed line), bcc (solid line), fcc (dashed-dotted line), and sh (dotted line) structures of Mg using the PAW pseudopotential and the fcc (stars) structure using the FP-LAPW method at 0 K. The open triangles, solid squares, and open squares are the measured data for hcp from refs 11, 32, and 12. The crosses are the experimental data for bcc from ref 11.

with only one atom in the primitive unit cell) is discovered with a ) 1.89 Å and c/a ) 0.90. We found that the c/a ratio in the sh structure is only weakly dependent on pressure (e.g., at 1000 GPa, c/a ) 0.91). It is noteworthy that the structural prediction performed at 1000 GPa also found the same sh structure, indicating that the stable field of this phase is at least up to 1000 GPa. It is noted that the sc phase proposed by Ahuja et al.14 has clearly a higher enthalpy than that of the sh phase below 1040 GPa, above which the sc structure becomes favorable. The enthalpy curves relative to bcc as a function of pressure for the chosen structures are presented in Figure 1a. It is clearly seen that the calculation correctly predicts the hcp f bcc transition at 53 GPa, in excellent agreement with the experimental results of 50 ( 6 GPa.11 Here, bcc is stable up to 456 GPa, above which the fcc structure has a lower enthalpy. The bcc f fcc transition pressure (456 GPa) is quite different from earlier results (180 GPa by GPT or 790 GPa by LMTO).13 The large discrepancy is not clear. We point out here that proper construction of the pseudopotential is crucial to the calculation of transition pressure. We then have employed an FP-LAPW calculation through the WIEN2K code25 to double confirm our results. It turns out that the bcc f fcc transition pressure through the full potential calculations, as plotted in the inset of Figure 1a, is about 451 GPa, in good accordance with the result of our pseudopotential calculations. The fcc structure is stable up to 756 GPa, and then the sh structure takes over. The high-pressure sh structure is experimentally observed earlier in elemental P,26 Si,27,28 and Ge.29 Phonon calculations give a criterion for the structural stability. We thus calculated phonon dispersion curves for the fcc and sh phases using the supercell method. No imaginary phonon frequencies are found in the pressure ranges

High-Pressure Crystal Structures and Behavior of Mg

J. Phys. Chem. C, Vol. 114, No. 49, 2010 21747

Figure 2. Calculated band structure (left panel) and partial electronic DOS plots (right panel) for fcc (a) at 500 GPa and sh (b) at 800 GPa. ELF for fcc (c) at 500 GPa and sh (d) at 800 GPa.

of 456-756 GPa for the fcc structure and above 756 GPa for the sh structure in the whole Brillouin zone, suggesting that both of the two structures are dynamically stable. The bcc f fcc phase transition is known for alkali metals and can be explained by the Bain deformation mechanism.30,31 It is found that the sh lattice can be easily derived from the fcc phase by distortion of the R and γ angles from 60° to 90° synchronously. The fcc structure (R ) β ) γ ) 60°) evolves to the sh structure when R and γ increase to 90°. This transition can also be viewed as an antiparallel shuffle movement of adjacent (001) planes along the [110]fcc direction. Accompanying the transformation from fcc to sh, the coordination number decreases from 12 to 8. The sh structure is denser than fcc because of the obvious shorter Mg-Mg distance in the sh structure. The calculated equation of states (EOS) for the hcp, bcc, and sh phases are shown in Figure 1b. The EOS results for both hcp and bcc phases are in excellent agreement with the experimental data.11,12,32 The volume changes at bcc f fcc (0.7%) and fcc f sh (2.7%) transformations (though rather small) characterize first-order transitions. It is known that, at very high density, the neglect of the core (1s and 2s) potential might influence the calculation results. To examine the validity of the adopted PAW pseudopotential at such high pressures, the EOS of the fcc and sh structures are also calculated by the FP-LAPW method (Figure 1b). Excellent mutual agreement between the PAW and FP-LAPW calculations over the entire pressure range studied supports the current methods and the PAW pseudopotential adopted. To probe the electronic properties and chemical bonding features of the fcc and sh phases, we calculated the electronic band structure, the partial electronic density of states (DOS), and the valence electron localization function (ELF) at 500 and 800 GPa, as shown in Figure 2. It can be clearly seen that the fcc (Figure 2a) and sh (Figure 2b) structures are still metallic by evidence of the energy bands crossing over the Fermi level. The total electronic DOS at the Fermi level for the fcc and sh structures is 0.22 and 0.17 eV-1, respectively, both of which

are much smaller than that (0.42 eV-1) of the hcp phase at ambient pressure. This suggests that the metallicity of the fcc and sh structures are relatively weaker than that of hcp. This is apparently in contrast to the increased metallicity observed in groups III-VII elements. Note that the deviation from “ideal” metallic behavior has also been found in group I elements under pressure. Particularly, sodium becomes a transparent wide-gap insulator at 200 GPa.5 The most intriguing feature of the partial DOS (Figure 2a,b) is the growing of the d state. The weakening of metallicity is attributed to the drop of 3d bands in energy relative to the 3p bands and the increased p-d hybridization upon compression. This phenomenon is substantially different from that of the low-pressure hcp structure, in which the bands are basically free-electron like, predominantly s and p in character, and the s-p hybridization is important for the metallicity of Mg. We also calculated the electronic DOS for the neighboring Na and Al at different pressures. The DOSs for Na, Mg, and Al are similar at ambient pressure with all having a significant amount of d characters at very high pressures, but quite different from each other for p-d hybridizations. The p-d hybridization in Na is more pronounced than those in Al and Mg. As a consequence, at ultrahigh pressures, Na turns to an insulator,5 but Al and Mg remain metallic, though relatively weaker than their ambient states. Thus, it is suggested that the increased d states and their hybridization with p orbitals in Na, Mg, and Al under high pressure play a critical role in determining the stable structures and the electronic states. ELF was introduced by Becke and Edgecombe as a “simple measure of electron localization in atomic and molecular systems”.33 As can be seen in Figure 2, we find that the valence electrons in the fcc structure are largely localized in the interstitials (here, ELF values are close to 0.92), similar to the insulating fcc-Ca at 18 GPa.8 However, what makes fcc-Mg unique is that it remains metallic. For the sh phase, the marked accumulation occurs only in the open interstitial regions and the maxima (ELF ≈ 0.90) are found at the positions on the sites 2d (1/3, 2/3, 1/2). In other words, if one sets the stacking of hexagonal layers as A, the ELF maxima are exactly located

21748

J. Phys. Chem. C, Vol. 114, No. 49, 2010

Li et al. calculated. It turns out that the temperature contribution does not affect the phase transition order, but slightly changes the transition pressures. The calculated phase boundaries for both bcc f fcc and fcc f sh transitions are found to have positive Clapeyron slopes (dP/dT ) ∆S/∆V), indicating the decreased vibrational entropy across the phase transition. The calculated ∆S are found to be -0.16 and -0.17 eV at 1000 K for bcc f fcc and fcc f sh transitions, respectively. Conclusions

Figure 3. P-T phase diagram of Mg as a function of pressure and temperature. The black solid lines show our predicted hcp-bcc, bcc-fcc, and fcc-sc transition boundaries. The red dashed-dotted line is the experimental melting curve from ref 12. The blue short dashed and green dashed lines represent the hcp-bcc boundary calculation from refs 40 and 39, respectively. The open square is the experimental hcp-bcc transition point measured by ref 11.

between two A layers. As the number of ionic cores is half of the interstitial electron density maxima, sh-Mg is an analogue of AlB2-type (Ω phase) structure: the ionic cores form the Al sublattice and the interstitial density maxima the B sublattice. The Ω structure shows two nonequivalent crystallographic sites, corresponding to the Wyckoff positions 1a and 2d sites, respectively. The Ω phase is observed at high pressure in pure Zr, Ti, and Hf and in some transition-metal alloys, such as Ti-V and Zr-Nb alloys.34,35 The phenomenon of increased localization of valence electrons as they are repelled from the ionic core regions into interstitial regions is also found in alkali and other alkaline-earth metal elements, such as Li,7,9 Na,5,36 K,6 and Ca.8 The fcc and sh structures can be described as electrides formed by Mg ionic cores and localized interstitial electrons, which we refer to as pseudoanions. The nearest Mg-Mg distances for the fcc and sh structures are 2.094 and 1.702 Å at 500 and 800 GPa, respectively, which implies strong valencevalence overlap (the 3s orbital radii in Mg are 1.279 Å) between neighboring Mg atoms. This may be the origin for valence electrons repelled into the lattice interstices. The valence electron localization of sh-Mg is not as strong as hp4-Na (ELF ≈ 1.0).5 The repulsion by core electrons in hP4-Na is forceful because of the strong core-valence overlap and even significant core-core overlap between neighboring Na atoms. This could be the reason why hP4-Na is an insulator and sh-Mg is a metal, although they adopt the similar hexagonal structure. We have carefully calculated the complete phase transition diagrams of Mg at high pressures and temperatures within the quasi-harmonic approximation,37,38 as shown in Figure 3. The calculated hcp f bcc transition pressure at room temperature is 49 GPa, in excellent agreement with the experimental value of 50 GPa.11 It can be clearly seen that our phase boundary results are in excellent agreement with experiment11 and earlier calculations39,40 below 90 GPa. We present here the higherpressure phase diagrams of Mg for the first time. Our results suggest that the transition pressures of bcc f fcc and fcc f sh increase slightly with increasing temperature. For example, the bcc f fcc transition pressure slightly increases from 461 GPa at 300 K to 478 GPa at 2000 K, a maximum temperature

In summary, we have extensively explored the high-pressure crystal structures of Mg by the PSO algorithm. In agreement with experiments, we reproduced the hcp and bcc structures at lower pressures. The predicted transition pressures, EOS, and temperature-dependent phase diagrams of the hcp and bcc structures are in excellent agreement with the experiments. A high-pressure fcc phase has been uncovered to be stable from 456 to 756 GPa. We also report an sh phase to be stable above 756 GPa. The fcc and sh structures are both metallic. The ELF indicated that the valence electrons of fcc and sh are mostly localized in the interstitial sites, similar to alkali and other alkaline-earth elements, such as Na, Li, K, and Ca. The predicted temperature-pressure phase diagram of Mg to 850 GPa is of fundamental interest and particularly important for improving our understanding of the high-pressure scale. Acknowledgment. We thank the National Natural Science Foundation of China (NSFC) under Grant Nos. 10874054, 11025418, and 91022029; the NSFC awarded Research Fellowship for International Young Scientists under Grant No. 10910263; and the China 973 Program under Grant No. 2005CB724400 for financial support. References and Notes (1) Donohue, J. The Structures of the Elements; Wiley: New York, 1974. (2) Nelmes, R. J.; Allan, D. R.; McMahon, M. I.; Belmonte, S. A. Phys. ReV. Lett. 1999, 83, 4081. (3) McMahon, M. I.; Bovornratanaraks, T.; Allan, D. R.; Belmonte, S. A.; Nelmes, R. J. Phys. ReV. B 2000, 61, 3135. (4) Kenichi, T. Phys. ReV. B 1994, 50, 16238. (5) Ma, Y.; Eremets, M.; Oganov, A. R.; Xie, Y.; Trojan, I.; Medvedev, S.; Lyakhov, A. O.; Valle, M.; Prakapenka, V. Nature 2009, 458, 182. (6) Marque´s, M.; Ackland, G. J.; Lundegaard, L. F.; Stinton, G.; Nelmes, R. J.; McMahon, M. I.; Contreras-Garcı´a, J. Phys. ReV. Lett. 2009, 103, 115501. (7) Ma, Y.; Oganov, A. R.; Xie, Y. Phys. ReV. B 2008, 78, 014102. (8) Oganov, A. R.; Ma, Y.; Xu, Y.; Errea, I.; Bergara, A.; Lyakhov, A. O. Proc. Natl. Acad. Sci. U.S.A. 2010, 107, 7646. (9) Neaton, J. B.; Ashcroft, N. W. Nature 1999, 400, 141. (10) Young, D. A. Phase Diagram of the Elements; University of California Press: Berkeley, CA, 1990. (11) Olijnyk, H.; Holzapfel, W. B. Phys. ReV. B 1985, 31, 4682. (12) Errandonea, D.; Meng, Y.; Husermann, D.; Uchida, T. J. Phys.: Condens. Matter 2003, 15, 1277. (13) McMahan, A. K.; Moriarty, J. A. Phys. ReV. B 1983, 27, 3235. (14) Ahuja, R.; Eriksson, O.; Wills, J. M.; Johansson, B. Phys. ReV. Lett. 1995, 75, 3473. (15) Umemoto, K.; Wentzxovitch, R. M.; Allen, P. B. Science 2006, 311, 983. (16) Wang, Y.; Lv, J.; Zhu, L.; Ma, Y. Phys. ReV. B 2010, 82, 094116. (17) Kennedy, J.; Eberhart, R. Particle Swarm Optimization; IEEE: Piscataway, NJ, 1995. (18) Eberhart, R.; Kennedy, J. A New Optimizer Using Particle Swarm Theory; IEEE: New York, 1995. (19) Yoshida, H.; Kawata, K.; Fukuyama, Y.; Takayama, S.; Nakanishi, Y. IEEE Trans. Power Syst. 2000, 15, 1232. (20) Parsopoulos, K. E.; Vrahatis, M. N. IEEE Trans. EVol. Comput. 2004, 8, 211. (21) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. ReV. Lett. 1996, 77, 3865. (22) Kresse, G.; Furthmuller, J. Phys. ReV. B 1996, 54, 11169.

High-Pressure Crystal Structures and Behavior of Mg (23) Blochl, P. E. Phys. ReV. B 1994, 50, 17953. (24) Kresse, G.; Joubert, D. Phys. ReV. B 1999, 59, 1758. (25) Blaha, P.; Schwarz, K.; Madsen, G. K. H.; Kvasnicka, D.; Luitz, J. WIEN2k, An Augmented Plane WaVe + Local Orbitals Program for Calculating Crystal Properties; Karlheinz Schwarz, Techn. Universita¨t Wien: Wien, Austria, 2001. (26) Akahama, Y.; Kobayashi, M.; Kawamura, H. Phys. ReV. B 1999, 59, 8520. (27) Hu, J. Z.; Spain, I. L. Solid State Commun. 1984, 51, 263. (28) Hu, J. Z.; Merkle, L. D.; Menoni, C. S.; Spain, I. L. Phys. ReV. B 1986, 34, 4679. (29) Olijnyk, H.; Sikka, S. K.; Holzapfel, W. B. Phys. Lett. 1984, 103A, 137. (30) Dmitriev, V. P.; Gufan, Y. M.; Tolédano, P. Phys. ReV. B 1991, 44, 7248.

J. Phys. Chem. C, Vol. 114, No. 49, 2010 21749 (31) Xie, Y.; Ma, Y.; Cui, T.; Li, Y.; Qiu, J.; Zou, G. New J. Phys. 2008, 10, 063022. (32) Clendenen, G. L.; Drickamer, H. G. Phys. ReV. 1964, 135, A1643. (33) Becke, A. D.; Edgecombe, K. E. J. Chem. Phys. 1990, 92, 5397. (34) Jameson, J. C. Science 1963, 140, 72. (35) Sikka, S. K.; Vohra, Y. K.; Chidambaram, R. Prog. Mater. Sci. 1982, 27, 245. (36) Neaton, J. B.; Ashcroft, N. W. Phys. ReV. Lett. 2001, 86, 2830. (37) Ma, Y.; Tse, J. S. Solid State Commun. 2007, 143, 161. (38) Gao, G.; Oganov, A. R.; Ma, Y.; Wang, H.; Li, P.; Iitaka, T.; Zhou, G. J. Chem. Phys. 2010, 133, 144508. (39) Mehta, S.; Price, G. D.; Alfe, D. J. Chem. Phys. 2006, 125, 194507. (40) Moriarty, J. A.; Althoff, J. D. Phys. ReV. B 1995, 51, 5609.

JP108136R