Cite this Article
Lv Chao-Jia, Liu Lei, Gao Yang, Liu Hong, Yi Li, Zhuang Chun-Qiang, Li Ying, Du Jian-Guo. Structural, elastic, and vibrational properties of phase H: A first-principles simulation. Chinese Physics B, 2017, 26(6): 067401  
Permissions
Structural, elastic, and vibrational properties of phase H: A first-principles simulation
Lv Chao-Jia1, Liu Lei1, 2, †, Gao Yang3, Liu Hong1, Yi Li1, Zhuang Chun-Qiang4, Li Ying1, 2, Du Jian-Guo1
Key Laboratory of Earthquake Prediction, Institute of Earthquake Science, China Earthquake Administration, Beijing 100036, China
National Key Laboratory of Shock Wave and Detonation Physics, Mianyang 621000, China
Department of Mechanical Engineering, Texas Tech University, Lubbock 79409, USA
Institute of Microstructure and Properties of Advanced Materials, Beijing University of Technology, Beijing 100124, China

 

† Corresponding author. E-mail: liulei@cea-ies.ac.cn

Abstract

Phase H (MgSiO H , one of the dense hydrous magnesium silicates (DHMSs), is supposed to be vital to transporting water into the lower mantle. Here the crystal structure, elasticity and Raman vibrational properties of the two possible structures of phase H with Pm and symmetry under high pressures are evaluated by first-principles simulations. The cell parameters, elastic and Raman vibrational properties of the Pm symmetry become the same as the symmetry at GPa. The symmetrization of hydrogen bonds of the Pm symmetry at GPa results in this structural transformation from Pm to . Seismic wave velocities of phase H are calculated in a range from 0 GPa to 100 GPa and the results testify the existence and stability of phase H in the lower mantle. The azimuthal anisotropies for phase H are , ( symmetry) and , ( symmetry) at 0 GPa, and increase to , ( symmetry) and , (Pm symmetry) at 30 GPa. The maximum direction for phase H is [101] and the minimum direction is [110]. The anisotropic results of seismic wave velocities imply that phase H might be a source of seismic anisotropy in the lower mantle. Furthermore, Raman vibrational modes are analyzed to figure out the effect of symmetrization of hydrogen bonds on Raman vibrational pattern and the dependence of Raman spectrum on pressure. Our results may lead to an in-depth understanding of the stability of phase H in the mantle.

PACS: 74.25.nd;62.50.-p;63.20.dk
Keyword:phase H;elastic properties;Raman properties;first principles
1. Introduction

The presence of water in the lower mantle has various functions, such as dramatically lowering the melting temperature of mantle rocks, causing arc magmatism, reducing the magma viscosity and density, enhancing magma migration, influencing the elastic and rheological properties of the mantle, and broadening seismic discontinuities.[1] It is now accepted that water can percolate through the Earth’s interior by subducting slabs through hydrous minerals, such as serpentine.[2,3] The serpentine can transform into several dense hydrous magnesium silicates (DHMSs), depending on the pressure and temperature conditions.[4–6] Since the first discovery of 10-Å phase,[7] dense hydrous magnesium silicates including 3.65-Å phase, phase A, phase B, superhydrous phase B, phase D, phase E, etc. have received considerable attention because of their presumed role in the water resources of the Earth.[8–10]

Phase D (MgSi has been considered to be the most thermodynamically stable phase under lower mantle conditions.[11,12] It has been discovered that phase D transforms into the phase H at 44 GPa.[13] Phase H was first announced to be stable up to 52 GPa and had a crystal structure similar to that of CaCl -type SiO , of which half of the Si atoms replace Mg atoms and hydrogen atoms are located at the edges of the vacant octahedral sites (Fig. 1). Early theoretical studies predicted the crystal structure of phase H to be monoclinic with symmetry.[14] Since then, phase H has been discussed to evaluate its structural stability, elastic properties and the possibility of forming a solid solution with -AlOOH in the lower mantle condition.[14–16] Yet, due to the fact that different analytical methods are adopted, the results of previous studies are controversial. Possible structure of phase H suggested by Bindi et al.[15] was based on energy-dispersive x-ray diffraction while that by Nishi et al.[16] was based on in situ single-crystal x-ray diffraction. Besides, the quality of diffraction data was not accurate enough to determine the positions of hydrogen atoms. On the other hand, the work provided by first principles calculations could not take the thermal vibration effect into consideration.[17] As a result, the structure parameters and symmetry of phase H are discrepant in different studies.

Fig. 1. (color online) Structures of phase H with and Pm symmetry. The hydrogen atoms are not located in the middle of two oxygen atoms in structure with Pm symmetry (right), however, hydrogen atoms in structure (left) are located on the right side in the middle of couples of oxygen atoms.

In this study, in order to further understand the stability of phase H under high pressure, the crystal structure, elastic constants, seismic velocities and anisotropies of the two different possible symmetry phase H’s (Pm and symmetry) are calculated by using first-principles simulations. The symmetry is based on the fully optimized structure by Tsuchiya and Mookherjee[17] while Pm symmetry possesses hydrogen atoms drawn at the Wyckoff position, which is inspired by previous single-crystal x-ray diffraction research.[16] Raman vibrational properties are also analyzed to evaluate the changes of the crystal structure with pressure.

2. Methods
2.1. Computational methods

The first-principles calculations were performed using density functional perturbation theory (DFPT),[18] density functional theory (DFT),[19,20] and the plane wave pseudopotential method, as implemented in the CASTEP codes.[21] The generalized gradient approximation (GGA) with PBE parameterization was adopted to describe the exchange-correlation functional.[22] GGA functional was tested to provide a better performance than local density approximation (LDA) for simulations of hydrogen bonding.[23,24] Norm-conserving pseudopotentials[25] were employed to model the electron–ion interaction. The electronic wave functions were expanded in plane-waves with a kinetic energy cutoff of 1200 eV. The irreducible Brillouin zone of the phase H is sampled on a 5 5 Monkhorst–Pack mesh.[26] Both energy cutoff and k-point mesh are selected after a series of convergence tests for the precision of the calculation. The tests show that the difference among lattice parameters and total energy are restrained under 0.01% at higher accuracy. The total energy convergence criterion of 5 a.u. was used in the self-consistent field calculations.

The structures of phase H used in this study possess monoclinic structures with and Pm symmetries which were proved to be more thermodynamically stable (Fig. 1). These two adopted structures have ordered Mg and Si arrangements which are essential for the calculation of polarizability and Raman spectra. However, the structures of synthetic experimental samples belong to orthorhombic structures with symmetry Pnnm[15] or with .[16] The difference from our calculation was that the structures found in experiments have disordered sequences of Mg and Si atoms in which the occupancy of Mg and Si for each center of the octahedron were both 50%. All structure parameters and atomic coordinates were fully relaxed to a static configuration (0 K) by using BFGS geometry optimization algorithm.[27] The elastic constants were determined by stress–strain relationships.[28] The magnitude of all applied strains was 0.003 and the linear relationship was ensured to be enough for this strain range. The Raman properties were calculated using the linear response density functional perturbation theory (DFTP).[29] Spatial derivatives of the macroscopic polarization were calculated numerically along the eigenvectors of each Raman active phonon mode by calculating the polarization for each displacement through using the linear response formalism.[30] Once these derivatives are known, the Raman cross-section through an appropriate averaging space can be calculated. More information about Raman spectra by density functional theory was listed by Porezag and Pederson[31] and Refson et al.[32]

2.2. Benchmark calculations

To assess the performances of the DFT and DFPT total-energy approach in our study, test calculations of the cell parameters and elastic constants of phase H were performed. The differences in cell volumes, lengths of a, b, and c axes between the calculated results and corresponding experimental values[16] are 0.1%, 1.2%, 0.7%, and 0.4% respectively, or 2.7%, 1.7%, 0.6%, and 0.7% compared with the other experimental results[15] (Table 1). Comparing our calculated results with the calculated results of Tsuchiya,[14] the tiny differences can be ignored. Since no experimental results of elastic constants are available, the calculated values of elastic constants[17] are adopted for comparison. The differences in , , , , , , , and between this work and previous study are 2.2%, 5.4%, 2.9%, 6.6%, 3.2%, 10.4%, 2.8%, and 6.7% at 30 GPa, respectively. Notably, the pseudopotential of magnesium in former research[17] is generated by the method of von Barth and Car. In contrast, here we adopted norm-conserving pseudopotentials for all atoms, which may lead to the discrepancies in elastic constants. Such close agreement between the experimental and calculated values of lattice parameters and elastic constants demonstrates the validity of this computational project.

Table 1.

Benchmark calculations results.

.
3. Results
3.1. Crystal structure

Ever since the first-principles prediction of this new mineral was proposed,[14] phase H has been observed in the two separate experimental research studies.[15,16] However, due to the difficulties in detecting the positions of hydrogen atoms and pressure-induced rapid amorphization of the sample,[15,16] the structure of phase H was roughly judged to be an orthorhombic system with symmetry ( -AlOOH-type structure) and orthorhombic system with symmetry Pnnm (CaCl type structure), respectively. However, theoretical studies proposed different crystal structures: monoclinic structures with Pm and symmetry, which have been proved to be more thermodynamically stable than the experimental ones. The major difference between the theoretical and experimental structures is the arrangement of Mg and Si atoms. Both experimental structures possess disordered distributions of cations (Mg, Si), while theoretical structures possess ordered distributions of Mg and Si in octahedral sites. Furthermore, the symmetry possesses symmetric hydrogen bonds, while the Pm symmetry possesses asymmetric hydrogen bonds. The Pm and symmetries slightly deviate from experiential structures and were provided to be precise enough for the simulation of phase H.[17] Therefore, we select the two structures from theoretical studies as our investigation objects.[14]

Geometry optimization has been conducted on both Pm and structures from 0 GPa to 100 GPa. The detailed cell parameters of the two structures are listed in Table 2 and Supplement Tables A1, A2, and A3 in Appendix A. Figure 2 shows the pressure-lattice parameters of the two structures under high pressures. The lengths of a, b, and c axes decrease at almost the same ratio with increasing pressure, at the same time the angle decreases to 90 . Two-stage changes of lattice parameters of the two structures with increasing pressure are observed. The variation tendencies for lattice parameters of the two structures approach each other from 0 GPa to 30 GPa, above which the two structures almost keep the same lattice parameters. A possible explanation to this result is the symmetrization of hydrogen bonds. Since the first prediction of symmetrization of the hydrogen bonds in ice VII,[33] a large number of research studies have been conducted to understand the pressure dependence of hydrogen bonds’ symmetrization.[34–36] It is widely accepted that the bonded O–H distance shows relevance to the distance in many O– geometry systems. In a weak hydrogen bond system with phase H, where distance is more than 2.6 Å, the hydrogen bond length is almost fixed at 1 Å. Once the bond is under high pressure, the distance would decrease to around 2.4 Å and the O–H length could increase to the half of the distance, where the hydrogen is located in the middle of two oxygen atoms. The hydrogen bonds of phase H with Pm symmetry become symmetrized at pressures higher than 30 GPa (Fig. 2). When the distance reduces with pressure increasing from 0 GPa to 30 GPa, the length of O–H bond also increases rather than falls. The length of O–H is kept at half of the distance of under higher pressures, which means that the hydrogen atom is located in the middle of the pair of oxygen atoms. This result indicates that phase H would turn from Pm structure with asymmetric hydrogen bonds into the structure with symmetric hydrogen bonds at about 30 GPa. It is also the reason that the two structures with different initial lattice parameters can become the same at pressures higher than 30 GPa.

Table 2.

Structural properties of phase H.

.
Fig. 2. (color online) Lattice parameters of phase H each as a function of pressure: (a) angle and a, b, and c axes of phase H, (b) cell volume of phase H, (c) symmetrization of hydrogen bonds in phase H with Pm symmetry. The black dashed line represents the function of . At about 30 GPa, and reach the half of .

The volume–pressure relationships (Fig. 2) for the two structures of phase H are well described by the finite strain formulation of the third order Birch–Murnaghan.[37] Values of bulk moduli are 156.2 GPa and 177.9 GPa, values of the pressure derivative of the zero pressure bulk modulus ) are 4.40 and 4.18 and the values of initial unit-cell volume are 59.3 {\AA } and 58.4 Å for Pm and structures, respectively. The differences in initial unit-cell volume, and , between the results in this work (Pm) and previous calculation results ( )[17] are 2.9%, 5.9%, and 1.6%, respectively.

3.2. Elasticity

The pressure-dependent elastic constants for the two phase H are shown in Fig. 3 and Supplement Table A3. All elastic constants, except for and both of which maintain around 0 GPa, increase with pressure rising, which is different from the previous result.[17] We think that the discrepancy of pseudopotential selected for magnesium might result in the differences in lattice parameters. The differences in lattice parameters in turn result in the discrepancies in elastic constants. For instance, the values of lattice parameter a in this work and the previous work[17] have a difference of 0.40% under 30 GPa, and this leads to a difference of 0.58% in elastic constant . The principle elastic constants and , the shear elastic constant , and the off-diagonal elastic constants , , and for Pm symmetry are significantly smaller than those of the symmetry under low pressures. These constants also sharply rise from 25 GPa to 30 GPa, which is also the pressure range of symmetrization of hydrogen bonds. These results indicate that the elastic constants of Pm symmetry are obviously dependent on the strength of the hydrogen bond. Namely, the symmetrization of hydrogen bonds increases the stiffness of phase H in the direction perpendicular to the c axis direction. Similar anomalous phenomena have been reported in phase D and -AlOOH.[35,36] Above 30 GPa, elastic constants for the two structures become equal to each other.

Fig. 3. (color online) Pressure-dependent elastic constants for phase H.

Bulk (K) and shear (G) moduli are determined by the Voigt–Reuss–Hill average (Fig. 4). The bulk and shear moduli of phase H with Pm symmetry show anomalous increases at 25 GPa–30 GPa. The anomalous pressure-dependent bulk and shear moduli are also related to the symmetrization of hydrogen bonds. The bulk modulus of phase H with Pm symmetry is about 27% smaller than that of symmetry. Under higher pressures ( GPa) the bulk and shear moduli of these two structures are almost the same.

Fig. 4. (color online) Bulk and shear moduli (a) and seismic velocities (b) of phase H as a function of pressure. Grey lines and points in panel (a) stand for K and G from Ref. [17].

Sound wave velocities from seismic observation can reflect the composition and structure of the mantle. For crystalline material, the compressional wave velocity and shear wave velocity are related to the bulk and shear moduli via

where is the mass density determined by the equilibrium of crystal lattice from the first-principles optimization, K is the bulk modulus, and G is the shear modulus. The and of phase H with symmetry increase smoothly with increasing pressure. The and of Pm symmetry sharply increase from 25 GPa to 30 GPa (Fig. 4) compared with those of . The pressure dependent and of phase H with Pm and symmetry show similar tendencies to those of K and G. At pressures above 30 GPa, the values of and of Pm and symmetry are almost identical. The symmetrization of hydrogen bonds also generally takes place in plenty of minerals under high pressures, such as ice, -AlOOH and hydrous silica.[36,38] Therefore, the symmetrization of hydrogen bonds in dense hydrous silicate may contribute to the abnormally pressure-dependent and in the upper mantle. The symmetry shows higher density and velocity than the Pm symmetry under pressure lower than 30 GPa. The transition from Pm to may somewhat contribute to the leaps of density and seismic wave velocity in the transition layer of the mantle.

The values of density, and of phase D, H, and -AlOOH under high pressures are plotted in Fig. 5. Phase D shows the minimal density and velocities while -AlOOH possesses the maximal density and velocities in the hull depth range. It has been proved that phase D would transform into phase H at pressures ranging from about 40 GPa to 50 GPa, and phase H is inclined to form solid solutions with -AlOOH while -AlOOH and H show similar structures.[15,16,39] The transformation process from phase D to phase H and finally to solid solution of phase H and -AlOOH results in the increasing of density and velocities, which show the consistent trends with PREM and AK135 patterns. -AlOOH is believed to be stable under lower mantle condition. The forming of solid solution with phase H implicates that the water may have an opportunity to be carried to the lower region of the mantle and even the border of the mantle and core of the Earth. Ultimately, water should be released freely after the solid solution of -AlOOH and phase H have been heated up by the core, which may contribute to the phenomena of ultralow-velocity regions[40] and large low-shear-velocity provinces[41] detected in the underpart of the mantle.

Fig. 5. (color online) Densities and seismic wave velocities of phase H, phase D and -AlOOH. Data for phase D are cited from Ref. [49] while data for AlOOH are taken from Ref. [50].

Single crystal seismic velocities for phase H are also calculated by solving Cristoffel’s equation:[42] , where n, , V, and are the propagation direction, density, elastic wave velocity, and Kronecker delta, respectively. Longitudinally and two orthogonally polarized shear wave velocities of the Pm and structure at 0, 30, and 60 GPa are shown in Fig. 6 by using the MSAT codes.[43] At 0 GPa, both and of Pm structure (solid lines) are significantly smaller than those of structure (dashed lines) except for those in the [001] direction, since there are only octahedrons of Si and O but no hydrogen bonds parallel to [001] direction. With the symmetrization of hydrogen bonds at 30 GPa, the patterns of single crystal seismic velocities for the two structures become similar. As pressure increases, the difference in anisotropy of seismic velocity between the two structures decreases significantly. The anisotropic properties of phase H at 30 GPa and 60 GPa are almost the same, so the primary reason why the anisotropy changes is crystal structural rearrangement (symmetrization of hydrogen bonds). The maximal direction for phase H is [101] and the minimal direction is [110]. Azimuthal anisotropy for P and S waves is defined as , where and are the maximal and minimal velocity, respectively. The azimuthal anisotropies for phase H are %, % ( symmetry) and %, % (Pm symmetry) at 0 GPa, and increase to %, % ( symmetry) and %, % (Pm symmetry) at 30 GPa. Those results show that the phase H can be a significant source of anisotropy in the lower mantle.

Fig. 6. (color online) Single crystal seismic velocities of phase H under 0, 30, and 60 GPa.
3.3. Vibration

Vibrational spectra of crystalline solids can be interpreted by using the factor group analysis.[44] For phase H with symmetry, the irreducible representation of optical vibrations is . Since the structure is centro-symmetric, IR and Raman-active modes are mutually exclusive.[45] The Raman-active modes are . As for Pm structure, the optical vibration mode is and Raman-active mode is . The Raman-active modes for experimental results (Pnnm and P21nm) are 3 and 3 , which deviate from calculated ones since different symmetries are involved. Nevertheless, the variance is quite acceptable as the major difference between experimental and calculated structures is the angle which distorts within a range of 3 degrees.

Raman frequencies of the two phase H structures are calculated and plotted each as a function of pressure in Fig. 8. The number of active Raman modes for symmetry remains constant in a pressure range from 0 GPa to 100 GPa. In addition, the vibrational frequencies increase with the augment of pressure, which is normal behavior resulting from the increased repulsion between neighbor atoms.[46] The intensities of all Raman-active modes except for decrease under pressure (Supplement Tables A4, A5, and A6 in Appendix A). As for symmetry, all Raman modes (Fig. 7) are only related to the motion of oxygen atoms. No contributions from the motions of Mg, Si, and H atoms are detected. The and modes only have atomic vibrations parallel to the c axis, and the , , , and modes only have atomic vibrations perpendicular to the c axis. Raman intensities for vibrations parallel to c axis ( and are smaller than those of vibrations perpendicular to c axis ( , , , and modes) (Supplement Table A4 in Appendix A). The results suggest that the Raman vibration in phase H with symmetric hydrogen bonds is mainly controlled by displacements of oxygen atoms parallel to (001) plane.

Fig. 7. (color online) Vibrational patterns of the structure under 0 GPa. Yellow, green, red, and white balls donate Si, Mg, O, and H atoms, respectively. The direction and scale of arrow indicate vibrational orientation and relative vibrational intensity.
Fig. 8. (color online) Plots of Raman-active frequency of phase H versus pressure. Dasheed plots and lines refer to the Raman modes which are inactive under high pressures.

At 0 GPa, 21 Raman modes are observed in phase H with Pm symmetry. These modes can be classified as two groups. Group 1 includes 15 vibrational modes with frequencies ranging from 226 cm to 780 cm . The rest of these modes with frequencies in a range from 1009 cm to 2470 cm belong to group 2. Distortions and relative displacements of Si (or Mg) and O are observed in group 1 while only displacement of H, stretching and bending of hydrogen bonds are observed in group 2. The O–H stretching modes ( and have the maximal intensities, because of the long-range electrostatic interactions between O–H groups.[41]

The vibrational patterns of Pm symmetry can be divided into two stages: from 0 GPa to 30 GPa and from 30 GPa to 100 GPa (Fig. 8), according to the symmetrization of hydrogen bonds at about 30 GPa. All the vibrational modes belonging to group 2 (in a frequency range from 1009 cm to 2479 cm ) disappear at pressures higher than 30 GPa. The reason is that the reduction of distances and the symmetrization of hydrogen bonds restrict the ability for hydrogen atoms to motion.[47] Besides, 8 vibrational modes in group 1 which relate to the motions of Mg or Si atoms vanish before being pressurized to 30 GPa. Finally, in the Pm symmetry there are only 6 vibration modes ( , , , , , and ) that remain Raman-active, the number of which modes (i.e., 6) is exactly the number of vibration modes that are Raman active in the symmetry. Despite the obvious differences between Raman spectra at 0 GPa, the Raman frequencies and intensities for Pm and symmetry are quite alike at pressures higher than 30 GPa (Supplement Tables A4, A5, and A6 in Appendix A).

Six Raman modes: , , , , , and of Pm symmetry are taken into account for calculating the pressure derivatives ( in a pressure range from 0 GPa to 100 GPa because the remaining modes vanish under high pressures ( GPa) (Table 1). The pressure derivatives of all modes for phase H are positive, which implies the strengthening of the structure.[48] Unlike the Raman shift frequency of symmetry which is linearly dependent on pressure, the pressure derivatives of Pm symmetry can be classified as 2 groups. At low pressures (< 25 GPa), the pressure derivatives of Raman frequencies for Pm symmetry are larger than those for symmetry. Also, the derivatives for both structures become similar at high pressures (>30 GPa). A conspicuous jump of frequencies for Pm structure takes place at pressures between 25 GPa to 30 GPa, which corresponds to symmetrization of hydrogen bonds.

Table 3.

Pressure shifts of the Raman modes of phase H.

.
4. Conclusions and perspectives

The phase H (MgSiO , as an important candidate for transporting water into the deep Earth, is stable at pressures greater than the stability field of phase D, and separately discussed to evaluate its structure stability, elastic properties and the possibility of forming a solid solution with -AlOOH in the lower mantle condition.[14–16] In order to understand the stability of phase H under high pressure, first-principles simulations are carried out to investigate the structural, elastic, and vibrational properties of phase H with the two possible symmetric structures: Pm and . The results indicate that the two structures show similar structural parameters, elastic and vibrational properties at pressures higher than 30 GPa. The reason is the symmetrization of hydrogen bonds. The symmetrization of hydrogen bonds also generally takes place in plenty of minerals under high pressure, such as ice, -AlOOH and hydrous silica.[36,38]

The and of phase H with Pm symmetry sharply increase from 25 GPa to 30 GPa compared with those with symmetry, and when pressure is larger than 30 GPa, the values of and of Pm and symmetry are almost identical. In other words, the symmetrization of hydrogen bonds made the seismic velocities of phase H increase. Therefore, the symmetrization of hydrogen bonds in dense hydrous silicate may lead to the abnormal pressure dependence of and and the formation of a high-velocity layer in the mantle. The azimuthal anisotropies for phase H are %, % ( ) and %, % (Pm) at 0 GPa, and increase to %, % ( symmetry) and %, % (Pm symmetry) at 30 GPa. So phase H might be one of the sources of the velocity anomalies and anisotropy in the upper part of the lower mantle. Compared with phase D and -AlOOH, phase H shows small density but fast seismic velocity, which might contribute to the anomalous phenomenon of seismic velocity in the lower mantle.

Raman spectra of phase H show the effects of symmetrization of hydrogen bonds on vibrational frequency and intensity. At 0 GPa, 21 Raman modes are observed in phase H with Pm symmetry. As for symmetry, 6 Raman modes are observed at 0 GPa. However, only 6 Raman modes in each phase H with Pm-to- symmetry are detected under higher pressure ( GPa). When pressure is higher than 30 GPa, Raman frequencies of Pm symmetry are the same as those of symmetry, which indicates the transformation from Pm to symmetry as a result of symmetrization of hydrogen bonds. The values of pressure derivative of all Raman modes are positive. All these results conduce to the understanding of the stability of phase H in the mantle.

Reference
[1] Huang X Xu Y Karato S 2005 Nature 434 746
[2] Faccenda M Burlini L Gerya T V Mainprice D 2008 Nature 455 1097
[3] Komabayashi T Omori S Maruyama S 2005 Phys. Earth Planet. Int. 153 191
[4] Akimoto S Akaogi M 1980 Phys. Earth Planet. Int. 23 268
[5] Ohtani E Toma M Litasov K Kubo T Suzuki A 2001 Phys. Earth Planet. Int. 124 105
[6] Ghosh S Schmidt M W 2014 Geochim. Cosmochim. Acta 145 72
[7] Sclar C B Carrison L C Schwartz C M 1965 Trans. Am. Geophys. Union 46 184
[8] Ohtani E Mizobata H Yurimoto H 2000 Phys. Chem. Miner. 27 533
[9] Komabayashi T Omori S Maruyama S 2004 J. Geophys. Res.: Solid Earth 109 B03206
[10] Melekhova E Schmidt M W Ulmer P Pettke T 2007 Geochim. Cosmochim. Acta 71 3348
[11] Liu L 1987 Phys. Earth Planet. Int. 49 142
[12] Frost D J Fei Y 1998 J. Geophys. Res.: Solid Earth 103 7463
[13] Shieh S R Mao H Hemley R J Li C M 1998 Earth Planet. Sci. Lett. 159 13
[14] Tsuchiya J 2013 Geophys. Res. Lett. 40 4570
[15] Bindi L Nishi M Tsuchiya J Irifune T 2014 Am. Mineral. 99 1802
[16] Nishi M Irifune T Tsuchiya J Tange Y Nishihara Y Fujino K Higo Y 2014 Nat. Geosci. 7 224
[17] Tsuchiya J Mookherjee M 2015 Sci. Rep. 5 15534
[18] Refson K Tulip P R Clark S J 2006 Phys. Rev. 73 155114
[19] Hohenberg P Kohn W 1964 Phys. Rev. 136 B864
[20] Kohn W Sham L J 1965 Phys. Rev. 140 A1133
[21] Clark S J Segall M D Pickard C J Hasnip P J Probert M J Refson K Payne M C 2005 Kristallogr. 220 567
[22] Perdew J P Wang Y 1992 Phys. Rev. 45 13244
[23] Hamann D R 1997 Phys. Rev. 55 R10157
[24] Umemoto K Wentzcovitch R M Baroni S De G S 2004 Phys. Rev. Lett. 92 105502
[25] Troullier N Martins J L 1991 Phys. Rev. 43 1993
[26] Monkhorst H J Pack J D 1976 Phys. Rev. 13 5188
[27] Pfrommer B G Cote M Louie S G Cohen M L 1997 J. Comput. Phys. 131 233
[28] Karki B B Stixrude L Wentzcovitch R M 2001 Rev. Geophys. 39 507
[29] Baroni S De G S Dal C A Giannozzi P 2001 Rev. Mod. Phys. 73 515
[30] Gonze X 1997 Phys. Rev. 55 10337
[31] Porezag D Pederson M R 1996 Phys. Rev. 54 7830
[32] Refson K Tulip P R Clark S J 2006 Phys. Rev. 136 864
[33] Holzapfel W B 1972 J. Chem. Phys. 56 712
[34] Pruzan P Chervin J C Wolanin E Canny B Gauthier M Hanfland M 2003 J. Raman Spectrosc. 34 591
[35] Tsuchiya J Tsuchiya T Tsuneyuki S 2005 Am. Mineral. 90 44
[36] Sano-Furukawa A Komatsu K Vanpeteghem C B Ohtani E 2008 Am. Mineral. 93 1558
[37] Murnaghan F D 1944 Proc. Natl. Acad. Sci. 30 244
[38] Goncharov A F Struzhkin V V Somayazulu M S Hemley R J Mao H K 1996 Science 273 218
[39] Ohira I Ohtani E Sakai T Miyahara M Hirao N Ohishi Y Nishijima M 2014 Earth Planet. Sci. Lett. 401 12
[40] Wen L Helmberger D V 1998 Science 279 1701
[41] Richards M A Duncan R A Courtillot V E 1989 Science 246 103
[42] Musgrave M J P 1970 Crystal Acoustics San Fransisco Holden-Day
[43] Walker A M Wookey J 2012 Comput. Geosci. 49 81
[44] Fateley W G McDevitt N T Dollish F R 1971 Appl. Spectrosc. 25 155
[45] Delattre S Balan E Lazzeri M Blanchard M Guillaumet M Beyssac O Haussuhl E Winkler B Salje E K H Calas G 2012 Phys. Chem. Miner. 39 93
[46] Huang X L Li F F Huang Y P Wu G Li X Zhou Q Liu B B Cui T 2016 Chin. Phys. 25 037401
[47] Umemoto K Wentzcovitch R M 2005 Phys. Rev. 71 012102
[48] Hemley R J Jephcoat A P Mao H K Zha C S Finger L W Cox D E 1987 Nature 330 737
[49] Mainprice D Page Y L Rodgers J Jouanna P 2007 Earth Planet. Sci. Lett. 259 283
[50] Tsuchiya J Tsuchiya T 2009 Phys. Earth Planet. Int. 174 122