## New Research In

### Physical Sciences

### Social Sciences

#### Featured Portals

#### Articles by Topic

### Biological Sciences

#### Featured Portals

#### Articles by Topic

- Agricultural Sciences
- Anthropology
- Applied Biological Sciences
- Biochemistry
- Biophysics and Computational Biology
- Cell Biology
- Developmental Biology
- Ecology
- Environmental Sciences
- Evolution
- Genetics
- Immunology and Inflammation
- Medical Sciences
- Microbiology
- Neuroscience
- Pharmacology
- Physiology
- Plant Biology
- Population Biology
- Psychological and Cognitive Sciences
- Sustainability Science
- Systems Biology

# Quantum Monte Carlo computations of phase stability, equations of state, and elasticity of high-pressure silica

Edited by Russell J. Hemley, Carnegie Institution of Washington, Washington, DC, and approved April 6, 2010 (received for review October 20, 2009)

## Abstract

Silica (SiO_{2}) is an abundant component of the Earth whose crystalline polymorphs play key roles in its structure and dynamics. First principle density functional theory (DFT) methods have often been used to accurately predict properties of silicates, but fundamental failures occur. Such failures occur even in silica, the simplest silicate, and understanding pure silica is a prerequisite to understanding the rocky part of the Earth. Here, we study silica with quantum Monte Carlo (QMC), which until now was not computationally possible for such complex materials, and find that QMC overcomes the failures of DFT. QMC is a benchmark method that does not rely on density functionals but rather explicitly treats the electrons and their interactions via a stochastic solution of Schrödinger’s equation. Using ground-state QMC plus phonons within the quasiharmonic approximation of density functional perturbation theory, we obtain the thermal pressure and equations of state of silica phases up to Earth’s core–mantle boundary. Our results provide the best constrained equations of state and phase boundaries available for silica. QMC indicates a transition to the dense *α*-PbO_{2} structure above the core-insulating D” layer, but the absence of a seismic signature suggests the transition does not contribute significantly to global seismic discontinuities in the lower mantle. However, the transition could still provide seismic signals from deeply subducted oceanic crust. We also find an accurate shear elastic constant for stishovite and its geophysically important softening with pressure.

## Introduction

Silica is one of the most widely studied materials across the fields of materials science, physics, and geology. It plays important roles in many applications, including ceramics, electronics, and glass production. As the simplest of the silicates, silica is also one of the most ubiquitous geophysically important minerals. It can exist as a free phase in some portions of the Earth’s mantle. In order to better understand geophysical roles silica plays in Earth, much focus is placed on improving knowledge of fundamental silica properties. Studying structural and chemical properties (1) offers insight into the bonding and electronic structure of silica and provides a realistic testbed for theoretical method development. Furthermore, studies of free silica under compression (2–7) reveal a rich variety of structures and properties, which are prototypical for the behavior of Earth minerals from the surface through the crust and mantle. However, the abundance of free silica phases and their role in the structure and dynamics of deep Earth is still unknown.

Free silica phases may form in the Earth as part of subducted slabs (8) or due to chemical reactions with molten iron (9). Determination of the phase stability fields and thermodynamic equations of state are crucial to understand the role of silica in Earth. The ambient phase, quartz, is a fourfold coordinated, hexagonal structure with nine atoms in the primitive cell (2). Compression experiments reveal a number of denser phases. The mineral coesite, also fourfold coordinated, is stable from 2–7.5 GPa, but we do not study it due to its large, complex structure, which is a 24 atom monoclinic cell (3). Further compression transforms coesite to a much denser, sixfold coordinated phase called stishovite, stable up to pressures near 50 GPa. Stishovite has a tetragonal primitive cell with six atoms (4). In addition to the coesite-stishovite transition, quartz metastably transforms to stishovite at a slightly lower pressure of about 6 GPa. Near 50 GPa, stishovite undergoes a ferroelastic transition to a CaCl_{2}-structured polymorph via instability in an elastic shear constant (5, 10–16). This transformation is second order and displacive, where motion of oxygen atoms under stress reduces the symmetry from tetrahedral to orthorhombic. Experiments (6, 7, 17) and computations (18–20) have reported a further transition of the CaCl_{2}-structure to an *α*-PbO_{2}-structured polymorph at pressures near the base of the mantle.

The importance of silica as a prototype and potentially key member among lower mantle minerals has prompted a number of theoretical studies (10, 13, 15, 16, 18–23) to investigate high-pressure behavior of silica. Density functional theory (DFT) successfully predicts many qualitative features of the phase stability (18–22), structural (21), and elastic (10, 13, 15, 16, 23) properties of silica, but it fails in fundamental ways, such as in predicting the correct structure at ambient conditions and/or accurate elastic stiffness (21, 22). In this work, we use the quantum Monte Carlo (QMC) method (24, 25) to compute silica equations of state, phase stability, and elasticity, documenting improved accuracy and reliability over DFT. Our work significantly expands the scope of QMC by studying the complex phase transitions in minerals away from the cubic oxides (26, 27). Furthermore, our results have geophysical implications for the role of silica in the lower mantle. Though we find the CaCl_{2}-*α*-PbO_{2} transition is not associated with any global seismic discontinuity, such as D”, the transition should be detectable in deeply subducted oceanic crust.

### DFT Method and Drawbacks.

DFT is currently the standard model for computing materials properties and has been successful in thousands of studies for a wide range of materials. DFT is an exact theory, which states that ground-state properties of a material can be obtained from functionals of the charge density alone. In practice, exchange-correlation functionals describing many-body electron interactions must be approximated. A common choice is the local density approximation (LDA), which uses a functional of the charge density in the uniform electron gas, calibrated via QMC (28). Alternate exchange-correlation functionals, such as generalized gradient approximations (GGAs) and hybrids are also useful, but there is no useable functional that can provide exact results.

DFT functionals are often expected to only have problems with materials exhibiting exotic electronic structures, such as highly correlated materials. However, silica has simple, closed shell, covalent/ionic bonding (1) and should be ideal for DFT. LDA provides good results for properties of individual silica polymorphs, but it incorrectly predicts stishovite as the stable ground-state structure rather than quartz. The GGA improves the energy difference between quartz and stishovite (21, 22), giving the correct ground-state structure, and this discovery was one of the results that popularized the GGA. However, almost all other properties of silica, such as structures, equations of state, and elastic constants are often worse within the GGA (21, 22) and other alternate approximations. Indeed, a drawback of DFT is that there exists no method to estimate the size of functional bias or to know which functional will succeed at describing a given property.

### QMC Method.

Explicit treatment of many-body electron interactions with QMC has shown great promise in studies of the electron gas, atoms, molecules, and simple solids (24, 25). Here we address the question of whether QMC is ready for accurate computations of complex minerals by studying silica. In QMC the only essential inputs are the trial wave function (we use a determinant of DFT orbitals) and pseudopotentials for the core electrons. Two common QMC algorithms are variational Monte Carlo (VMC) and diffusion Monte Carlo (DMC). The VMC method is efficient and provides an upper bound on the ground-state energy using the variational principle applied to a trial form of the many-body wave function. However, VMC usually cannot give the required chemical accuracy. DMC refines the many-body wave function statistically and projects out the ground-state solution consistent with the nodes of the trial wave function. All QMC computations reported here are DMC results unless stated otherwise. Computed results, such as the total energy, have a statistical uncertainty that decreases as , where *N* is the number electronic configurations sampled. We propagated statistical errors into thermodynamic quantities through a linearized Taylor expansion or Monte Carlo techniques.

## Results and Discussion

### Equations of State.

Fig. 1 shows the computed equations of state compared with experimental data for quartz (2, 29), stishovite/CaCl_{2} (30, 31), and *α*-PbO_{2} (6, 7). We compute thermal equations of state from the Helmholtz free energy (32), *F*(*V*,*T*) = *U*_{static}(*V*,*T*) - *TS*_{vib}(*V*,*T*). We use QMC to compute the internal energy of the static lattice, *U*_{static}, and DFT linear response in the quasiharmonic approximation provides the vibrational energy contribution, *TS*_{vib}, where *T* is temperature and *S*_{vib} is vibrational entropy. In general, there is also an electronic contribution to the entropy due to thermal excitations in materials with small band gaps or metals. Since silica is a large band gap insulator, electronic entropy may be ignored for temperature ranges we consider. We fit isotherms of the Helmholtz free energy using an analytic equation of state with the Vinet form (33). Pressure is determined from *P* = -(∂*F*/∂*V*)_{T}. The gray shading of the QMC curves indicates one standard deviation of statistical error from the Monte Carlo data. The QMC results generally agree well with diamond-anvil-cell measurements at room temperature, as do the corresponding DFT calculations using the GGA functional of Wu and Cohen (WC) (34).

We computed thermodynamic parameters from the thermal equations of state. Table 1 summarizes ambient computed and available experimentally measured values (6, 29–31, 35–38, 39) of the Helmholtz free energy, *F*_{0} (Ha/SiO_{2}), the Helmholtz free energy difference relative to quartz, Δ*F* (Ha/SiO_{2}), volume, *V*_{0} (Bohr^{3}/SiO_{2}), bulk modulus, *K*_{0}, pressure derivative of the bulk modulus, , thermal expansivity, α (K^{-1}), heat capacity, *C*_{p}/*R*, Grüneisen ratio, γ, volume dependence of the Grüneisen ratio, *q*, and the Anderson–Grüneisen parameter, *δ*_{T}, for the quartz, stishovite, and *α*-PbO_{2} phases of silica. QMC generally offers excellent agreement with experiment for each of these parameters.

### Enthalpy Differences.

A phase transition occurs when the difference in free energy (or enthalpy at *T* = 0) changes sign, as shown in Fig. 2. Statistical uncertainty in the energy differences determines how well phase boundaries are constrained. The one-sigma statistical error on the QMC enthalpy difference is 0.5 GPa for the quartz-stishovite transition and 8 GPa for the CaCl_{2}-*α*-PbO_{2} transition. The error on the latter is larger because the scale of the enthalpy difference between the quartz and stishovite phases is about a factor of 10 larger than for CaCl_{2} and *α*-PbO_{2}. In both transitions, variation in the DFT result with functional approximation is large. For the metastable quartz-stishovite transition, LDA incorrectly predicts stishovite to be the stable ground state, WC underestimates the quartz-stishovite transition pressure by 4 GPa, and the GGA of Perdew, Burke, and Ernzerhof (PBE) (40) matches the QMC result. For the CaCl_{2}-*α*-PbO_{2} transition, the same three DFT approximations lie within the statistical uncertainty of QMC. The variability of the present calculations is less than different experimental determinations of this transition (Fig. 3). The experimental variability may be due to the difficulty in demonstrating rigorous phase transition reversals as well as pressure and temperature gradients and uncertainties in state-of-the-art experiments.

### Improvement of DFT.

One of the goals of this work is to find what is needed to improve DFT functionals. Previous analysis of LDA and GGA charge densities and densities of states showed no changes in the predicted bonding of quartz and stishovite. However, changes were apparent in the atomic regions, akin to slightly different predictions of effective ionic size (1). Consistent with this analysis, we find that, for each phase, QMC energies and pressures are shifted by roughly a constant relative to DFT (WC) at every volume. For instance, QMC gives a lower energy (pressure) of 0.059 Ha/SiO_{2} (1.5 GPa) for quartz, 0.047 Ha/SiO_{2} (3.8 GPa) for stishovite, and 0.044 Ha/SiO_{2} (6.2 GPa) for *α*-PbO_{2} for any given cell volume. The constant shift suggests one could correct DFT computations for silica and potentially other systems by simply doing a couple QMC calculations for each phase.

### Phase Boundaries.

Fig. 3*A* compares QMC and DFT predictions with measurements for the quartz-stishovite phase boundary. The QMC boundary agrees well with thermodynamic modelling of shock data (41, 42) and thermocalorimetry measurements (39, 43), while the WC boundary is about 4 GPa too low in pressure. The melting curve shown is from a classical model (44), which agrees well with available experiments collected in the reference. The triple point seen in the melting curve is for the coesite-stishovite transition, and not the metastable quartz-stishovite transition that our computed boundaries represent. The geotherm (45) is shown for reference.

We show similar QMC and DFT calculations compared with experiments for the CaCl_{2}-*α*-PbO_{2} phase boundary in Fig. 3*B*. The WC boundary lies within the QMC statistical error. Previous DFT work also shows that the LDA boundary lies near the upper range of our QMC boundary and PBE produces a boundary 10 GPa higher than the LDA (19). The two diamond-anvil-cell experiments (6, 7) shown constrain the transition to lie between 65 and 120 GPa near the mantle geotherm (2500 K), while QMC constrains the transition to 105 (8) GPa. The boundary inferred from shock data (17) agrees well with QMC and WC. The boundary slope measured by Dubrovinsky et al. (6) is negative, which is in contrast to the positive slope inferred by Akins et al. (17). QMC and WC as well as previous DFT studies (19, 20) predict a positive slope.

### Geophysical Significance.

The QMC CaCl_{2}-*α*-PbO_{2} boundary indicates that the transition to *α*-PbO_{2}, within a two-sigma confidence interval, occurs in the depth range of 2,000–2,650 km (86–122 GPa) and in the temperature range of 2,300–2,600 K in the lower mantle. This places the transition 50–650 km above the D” layer, a thin boundary surrounding Earth’s core ranging from a depth of ∼2,700 to 2,900 km (46). The DFT boundaries all lie within the QMC two-sigma confidence interval, with PBE placing the transition most near the D” layer. Free silica in D”, such as in deeply subducted oceanic crust or mantle–core reaction zones, would have the *α*-PbO_{2} structure. However, based on our results, the absence of a global seismic anomaly above D” suggests that there is little or no free silica in the bulk of the lower mantle. The *α*-PbO_{2} phase is expected to remain the stable silica phase up to the core–mantle boundary. With respect to postperovskite, MgSiO_{3}, measurements (47, 48) at 120 GPa in D”, QMC predicts *α*-PbO_{2} has 1–2% lower density and 6–7% larger bulk sound velocity, which may provide enough contrast to be seen seismically if present in appreciable amounts.

### Shear Softening.

Most of our information about the deep Earth comes from seismology, and we need elastic constants to compute sound velocities in the Earth. Much work has been done using DFT to compute and predict elastic constants for minerals in the Earth (23), but there is much uncertainty in the predicted elastic constants because different density functionals give significantly different values. Here, we test the feasibility of using QMC for computing the softening of a shear elastic constant, *c*_{11}-*c*_{12}, in stishovite, which drives the ferroelastic transition to CaCl_{2} (5, 10–16). We determine the elastic constant by computing the total energy as a function of volume-conserving strains of up to 3% strain. Accurate computation of the QMC total energies on a small strain scale is very computationally expensive, requiring roughly 100–1,000 times more CPU time than a corresponding DFT calculation.

Fig. 4 shows our computations for *c*_{11}-*c*_{12}, which took over 3 million CPU hours. For a well-optimized trial wave function, VMC often comes close to matching the results of DMC. We therefore chose to compute *c*_{11}-*c*_{12} at several pressures using VMC and checked only the endpoints with the more accurate DMC method. We also show the result of our WC and previous LDA computations (10). Both QMC and DFT results correctly describe the softening of *c*_{11}-*c*_{12}, indicating the zero temperature transition to CaCl_{2} near 50 GPa. Radial X-ray diffraction data (11) lies lower than calculated results. However, discrepancies can arise in the experimental analysis depending on the strain model used. Recent Brillouin scattering data (14) agrees well with DMC.

## Conclusions

We have presented QMC computations of silica equations of state, phase stability, and elasticity. We find the CaCl_{2}-*α*-PbO_{2} transition is not associated with the global D” discontinuity, indicating there is not significant free silica in the bulk lower mantle. However, the transition should be detectable in deeply subducted oceanic crust. Additionally, we document the improved accuracy and reliability of QMC relative to DFT. DFT currently remains the method of choice for computing material properties because of its computational efficiency, but we have shown that QMC is feasible for computing thermodynamic and elastic properties of complex minerals. DFT is generally successful but does display failures independent of the complexity of the electronic structure and sometimes shows strong dependence on functional choice. With the current levels of computational demand and resources, one can use QMC to spot-check important DFT results to add confidence at extreme conditions or provide insight into improving the quality of density functionals. In any case, we expect QMC to become increasingly important and common as next generation computers appear and have a great impact on computational materials science.

## Methods

### Pseudopotentials.

In order to improve computational efficiency, pseudopotentials replace core electrons of the atoms with an effective potential. We constructed optimized nonlocal, norm-conserving pseudopotentials with the OPIUM (49) code. For DFT calculations, we generated each pseudopotential using the corresponding exchange-correlation functional (LDA, PBE, WC). QMC calculations used pseudopotentials generated with the WC functional. In all cases, the silicon potential has a Ne core with equivalent 3s, 3p, and 3d cutoffs of 1.80 a.u. The oxygen potential has a He core with 2s, 2p, and 3d cutoffs of 1.45, 1.55, and 1.40 a.u., respectively.

### DFT Computations.

We performed static DFT computations with the ABINIT package (50). All calculations used a plane-wave energy cutoff of 100 Ha and Monkhorst-Pack k-point meshes such that total energy converged to tenths of milli-Ha/SiO_{2} accuracy. We used *k*-point mesh sizes of 4 × 4 × 4, 4 × 4 × 6, and 4 × 4 × 4 for quartz, stishovite/CaCl_{2}, and *α*-PbO_{2}, respectively, and all were shifted by (0.5, 0.5, 0.5) from the origin. For each phase, we computed total energies for six or seven cell volumes ranging from roughly 10% expansion to 30% compression about the equilibrium volume. Constant-volume optimization of cell geometries relaxed forces on the atoms to less than 10^{-4} Ha/Bohr. We fit the Vinet equation of state to the corresponding energy points as a function of volume.

### Phonon Computations.

We computed phonon free energies with the ABINIT package. We modeled lattice dynamics using the linear response method (51) within the quasiharmonic approximation (32). We computed phonon free energies over a large range of temperatures for each cell volume and added them to ground-state energies in order to obtain equations of state at various temperatures. Phonon energies were computed up to the melting temperature in steps of 5 K. Converged calculations used a plane-wave cutoff energy of 40 Hartree with matching 4 × 4 × 4 *q*-point and *k*-point meshes.

### QMC Computations.

We performed QMC computations using the CASINO Code (52). The QMC calculations are composed of three major steps: (*i*) DFT calculation producing a relaxed crystal geometry and single particle orbitals, (*ii*) construction of a trial wave function and optimization within VMC, and (*iii*) a DMC calculation to determine the ground-state wave function accurately.

We construct a trial QMC wave function by multiplying the determinant of single particle DFT orbitals with a Jastrow correlation factor (24, 25). As a check for dependence on DFT functional choice, we compared QMC total energies for stishovite using various functionals for the orbitals and found the energies were equivalent within one-sigma statistical error (tenths of milli-Ha). The Jastrow describes various correlations by including two-body (electron–electron, electron–nuclear), three body (electron–electron–nuclear), and plane-wave expansion terms, with a total of 44 parameters. We optimized parameters in the Jastrow by minimizing the variance of the VMC total energy over several hundred thousand electron configurations. Finally, we performed DMC calculations using well-optimized trial wave functions. A typical DMC calculation used 4000 electron configurations, a time step of 0.003 Ha^{-1}, and several thousand Monte Carlo steps.

In order for QMC calculations of solids to be feasible, a number of approximations (24, 25) are usually implemented that can be classified into controlled and uncontrolled approximations. The controlled approximations include statistical Monte Carlo error, numerical grid interpolation (∼5 points/*Å*) of the DFT orbitals (53), finite system size, the number of configurations used, and the DMC time step. All of the controlled approximations are converged to within at least 1 milli-Ha/SiO_{2}. Finite size errors are reduced by using a model periodic Coulomb Hamiltonian (54) while simulating cell sizes up to 72 atoms (2 × 2 × 2) for quartz, 162 atoms (3 × 3 × 3) for stishovite/CaCl_{2}, and 96 atoms (2 × 2 × 2) for *α*-PbO_{2}. We correct additional finite size errors due to insufficient *k*-point sampling based on a converged *k*-point DFT calculation. The uncontrolled approximations include the pseudopotential, nonlocal evaluation of the pseudopotential, and the fixed node approximation. These errors are difficult to estimate, but we used the scheme of Casula (55) to minimize the nonlocal pseudopotential error and some evidence suggests the fixed node error may be small (56, 57).

## Acknowledgments

We acknowledge Ken Esler, Richard Hennig, and Neil Drummond for helpful discussions. We are grateful to Richard Hennig for QMC tests showing independence of DFT functional choice. This research is supported by the National Science Foundation (NSF) (EAR-0530282, EAR-0310139) and the Deparment of Energy (DOE) (DE-FG02-99ER45795). A portion of the computational resources were provided through the Breakthrough Science Program of the Computational and Information Systems Laboratory of the National Center for Atmospheric Research, which is supported by the NSF. This research also used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U. S. DOE under Contract No. DE-AC02-05CH11231. Additional computational resources were provided by the TeraGrid, National Center for Supercomputing Applications, Ohio Supercomputer Center, and Computational Center for Nanotechnology Innovations.

## Footnotes

^{1}To whom correspondence should be addressed. E-mail: driver{at}mps.ohio-state.edu.Author contributions: R.E.C. and B.M. designed research; K.P.D., R.E.C., Z.W., and B.M. performed research; K.P.D., R.E.C., and B.M. analyzed data; P.L.R., M.D.T., and R.J.N. contributed analytic tools; and K.P.D., R.E.C., and J.W.W. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.0912130107/-/DCSupplemental.

## References

- ↵
- Wright K,
- Catlow CRA

- Cohen RE

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Knittle E,
- Jeanloz R

- ↵
- Syono Y,
- Manghnani MH

- Cohen RE

_{2}and geophysical implications, eds Syono Y, Manghnani MH (Terra Scientific, Tokyo; American Geophysical Union, Washington, DC), pp 425–431. - ↵
- ↵
- Brazhkin VV,
- et al.

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Demuth Th,
- Jeanvoine Y,
- Hafner J,
- Ángyán JG

- ↵
- ↵
- ↵
- ↵
- Needs RJ,
- Towler MD,
- Drummond ND,
- López Ríos P

- ↵
- ↵
- ↵
- ↵
- Hazen RM,
- Finger LW,
- Hemley RJ,
- Mao HK

- ↵
- Andrault D,
- Fiquet G,
- Guyot F,
- Hanfland M

- ↵
- Andrault D,
- Angel RJ,
- Mosenfelder JL,
- Le Bihan T

- ↵
- Anderson OL

- ↵
- ↵
- ↵
- Carmichael RS

- Sumino Y,
- Anderson OL

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Sidorin I,
- Gurnis M,
- Helmberger DV

- ↵
- Murakami M,
- Hirose K,
- Kawamura K,
- Sata N,
- Ohishi Y

_{3}. Science 304:855–858. - ↵
- Hutko AR,
- Lay T,
- Revenaugh J,
- Garnero EJ

- ↵For OPIUM pseudopotential generation programs, see http://opium.sourceforge.net .
- ↵
- ↵
- ↵
- Needs RJ,
- Towler MD,
- Drummond ND,
- López Ríos P

- ↵
- Alfè D,
- Gillan MJ

- ↵
- ↵
- Casula M

- ↵
- ↵