Ab initio calculations for the $\beta$-tin diamond transition in silicon: Comparing theories with experiments

We investigate the pressure-induced metal-insulator transition from diamond to $\beta$-tin in bulk silicon, using quantum Monte Carlo (QMC) and density functional theory (DFT) approaches. We show that it is possible to efficiently describe many-body effects, using a variational wave function with an optimized Jastrow factor and a Slater determinant. Variational results are obtained with a small computational cost and are further improved by performing diffusion Monte Carlo calculations and an explicit optimization of molecular orbitals in the determinant. Finite temperature corrections and zero-point motion effects are included by calculating phonon dispersions in both phases at the DFT level. Our results indicate that the theoretical QMC (DFT) transition pressure is significantly larger (smaller) than the accepted experimental value. We discuss the limitation of DFT approaches due to the choice of the exchange and correlation functionals and the difficulty of determining consistent pseudopotentials within the QMC framework, a limitation that may significantly affect the accuracy of the technique.

Figure 1

Accuracy of the DFT energy as a function of the basis set extension $n(=n_0)$ defined in the text. The number of Gaussians in each angular channel $l$ is given by $n_l=n_0-2l$, as explained in the text. In the plot, we indicate the value of $n$ when the corresponding angular channel has been switched on in the basis set expansion. The red open squares refer to the metallic $\beta$-tin phase at volume 13.081 Å${}^3/\mathrm{Si}$ and $c/a=0.54$, while the blue triangles to the diamond insulating phase at volume 20.036 Å${}^3$/Si. The green circles are the energy difference between the two phases.

Figure 2

Finite-size extrapolation for the $\beta$-tin ($V=15$ Å${}^3$/Si, $c/a=0.55$) and the diamond ($V=19.949$ Å${}^3$/Si) phases in the upper and lower panel, respectively. The data points correspond to the values reported in Table , obtained with the Trail-Needs pseudopotentials. In the plot the energy per atom is shown as a function of $1/N_{\mathrm{atoms}}$, where $N_{\mathrm{atoms}}$ is the number of Si in the crystal supercell. Note the slopes for a given boundary condition (either PBC or APBC) in the two cases differ by a factor of 10, making the extrapolation much harder for the metallic phase.

Figure 3

EOS obtained by DFT LDA calculations with PWSCF and LAPW methods. The PPs in the PW formalism have been generated by the Troullier-Martins scheme with the cutoff radius $r_c$ reported in the figure key. The zero of the energy has been chosen such that all $\beta$-tin EOS curves collapse at the volume $V=16{\AA }^3$/Si. This choice helps the comparison among the EOS curves. Note that the PWSCF calculations done with the PAW pseudopotential are on top of the all-electron LAPW points.

Figure 4

Effect of harmonic quantum corrections on the transition pressure in silicon as a function of temperature.

Figure 5

The corrected equation of state once the zero-point energy and the temperature effects are added (estimated at the DFT-PBE level). The reference is the VMC calculation reported in Fig. 7, and here the zero of energy is taken as the VMC value at $V=13.08$ Å${}^3$/Si.

Figure 6

Energy versus $c/a$ in the $\beta$-tin phase at fixed volume $15$ Å${}^3$/Si for various sizes. FS corrections from the KZK extrapolation are included and energies are averaged over the two most symmetric points ($\Gamma$ and $M$) of the Brillouin zone. In the largest size the Jastrow factor was parametrized at large distance according to Eq. (2) and the error of the restricted variational freedom ($\simeq -0.0082$ eV/atom) was estimated with the same calculation for 64 Si atoms.

Figure 7

Equation of state and the metal-insulator transition pressure in bulk silicon obtained by the VMC technique after optimizing a Jastrow factor in a localized basis set containing $2s$ and $2p$ Gaussian orbitals per Si. The energies are corrected using the KZK correction scheme.[23] In this way finite-size effects of the transition pressure appears to be small. The lattice value $c/a=0.55$ is used for the $\beta$-tin phase. The data plotted here are not yet corrected for the quantum and temperature effects. See Table  for the final pressure which includes zero-point motion and thermal contributions.

Figure 8

The same as in Fig. 7 but for the LRDMC calculations.