Excited state calculations in solids by auxiliary-field quantum Monte Carlo
We present an approach for ab initio many-body calculations of excited states in solids. Using auxiliary-field quantum Monte Carlo, we introduce an orthogonalization constraint to prevent collapse of the stochastic Slater determinants in the imaginary-time propagation. Trial wave functions from density-functional calculations are used for the constraints. Detailed band structures calculated for standard semiconductors are in good agreement with and experimental results. For the challenging ZnO wurtzite structure, we obtain a fundamental band gap of eV, consistent with experiments.
pacs:71.15.Qe, 02.70.Ss, 71.20.-b, 71.10.-w
The accurate calculation of excited states in extended systems is a leading challenge in modern electronic structure theory. Density functional theory (DFT), which usually gives good results for ground state properties of materials without strong electron correlation, suffers from the well-known band-gap problem even for simple semiconductors Hohenberg and Kohn (1964); Kohn and Sham (1965); Perdew and Levy (1983). The use of hybrid functionals, which incorporate a portion of exact exchange from Hartree-Fock (HF) theory, has led to significant improvements for semiconductors with small to medium gaps, but not for large gap systems Becke (1993). Time-dependent DFT and many-body perturbative approaches have shown considerable promise Onida and Rubio (2002). The method is perhaps the simplest and most widely used of the latter and has led to dramatically improved results, largely for simple bonded materials Hedin (1965). It is less successful in materials that are more strongly correlated. In wurtzite structure ZnO, for example, the accuracy of quasi-particle excitation energies is still a matter of some controversy, involving many factors such as the choice of pseudopotential, approximate exchange-correlation functional, and choice of plasmon-pole model Shih et al. (2010); Stankovski et al. (2011); Partoens et al. (2010). A general approach which allows accurate calculations of electronic excitations across a wide variety of solids is thus very much in need.
Quantum Monte Carlo (QMC) Foulkes et al. (2001); Zhang and Krakauer (2003) is a non-perturbative, many-body computational method which is uniquely capable of scaling up to large system sizes and which in principle does not involve empirical parameters. Ground-state calculations with several flavors of QMC have seen significant recent development Foulkes et al. (2001); Zhang and Krakauer (2003); Booth et al. (2009). Although some excited states calculations have also been performed, for example, calculations in molecules Schautz et al. (2004); Purwanto et al. (2009a) and diffusion Monte Carlo (DMC) calculations of excitation energies of silicon and diamond at some high symmetry -points Williamson et al. (1998); Towler et al. (2000), the capability of QMC to treat excited states in general is much less developed. A major reason for this is the intrinsic difficulty of maintaining orthogonality with the lower-lying states when the targeted many-body excited state is being represented stochastically in an imaginary-time projection method.
The auxiliary-field quantum Monte Carlo (AFQMC) method Zhang and Krakauer (2003) offers a new framework for addressing this difficulty and doing excited state calculations in solids. The random walks in AFQMC take place in a space of Slater determinants. The single-particle orbitals in each Slater determinant are expressed explicitly in terms of a chosen one-particle basis. Thus, in addition to orthogonalizing the single-particle orbitals with each other (as is done in ground-state calculations), one can orthogonalize them with empty (hole) orbitals to remove contamination in the projection. For example, as an approximate constraint, the occupied orbitals can be orthogonalized with unoccupied orbitals defined by a trial wave function. Recent developments in ground state calculations have demonstrated that the AFQMC framework, with properly formulated sign or phase constraints using trial wave functions, allows sign-problem-free simulations with high accuracy across a wide range of both strongly correlated model Hamiltonians Chang and Zhang (2008) and real material systems Zhang and Krakauer (2003); Suewattana et al. (2007); Kwee et al. (2008); Al-Saidi et al. (2006).
In this paper, we show how the AFQMC approach can be formulated to accurately calculate excited states in extended systems. The ground-state method is augmented by an orthogonalization constraint with virtual orbitals to prevent the collapse to lower-lying states in the imaginary-time projection. Simple trial wave functions from HF or DFT are used for the phaseless and orthogonalization constraints. We illustrate the method by calculating the detailed band structures in diamond structure silicon and carbon, which are, respectively, small and large band gap semiconductors. We then apply the method to determine the bandgap in ZnO, which has drawn much recent attention and which presents a significant challenge, with high-lying strongly correlated -bands interacting with and semicore states.
In AFQMC, we target an excited state of the many-body Hamiltonian with the imaginary-time projection , where is a trial wave function of the excited state, and with the simulation time step and the Trotter step size. The wave function at any imaginary-time can be thought of as . Each Slater determinant in the sum is a random walker, and their distribution gives a statistical representation of the coefficients . Explicitly, has the form , where each is an orbital (expressed in terms of the chosen single-partle basis) that evolves with , and denotes the number of spin- electrons. (Below we will suppress the spin index where possible; the spin- component has a similar form.) Thus the AFQMC process resembles the propagation of a population of mean-field solutions in time-dependent external fields. Periodically in the propagation the orbitals within each random walker are orthogonalized with each other to ensure that the signal for fermionic antisymmetry is not lost in numerical noise. This step is also needed in ground-state calculations, and indeed even in mean-field calculations. The excited state energy can be calculated by the mixed estimator Purwanto et al. (2009a):
which converges to the exact result at . If belongs to an irreducible representation of the symmetry group of different from that of the ground state, Eq. (1) is exact for the lowest excited state of that symmetry. Otherwise, imaginary-time projection is considerably more challenging, since the propagation will tend to collapse to the ground state (or other lower-lying states) Purwanto et al. (2009a).
The AFQMC formalism, however, provides the ability to prevent the collapse by imposing an additional orthogonality constraint naturally, using virtual orbitals. For a concrete illustration, let us consider a targeted many-body state corresponding to the “single” excitation of replacing the valence orbital by a conduction orbital labeled by (thus ). Each random walker is still an -electron Slater determinant. For the purpose of orthogonalization, however, we will regard it as the extended ordered list . Any orbital denoted by ‘-’ is a virtual orbital whose only role is in the orthogonalization of the occupied orbitals. As an approximation, we will use the corresponding orbitals in , i.e., . The choice of is further discussed below. A Gram-Schmidt orthogonalization on this extended set ensures that the following orthogonality conditions are obeyed: i) ; ii) for ; iii) for . If the valence state is degenerate, its partners are not included in the constraint, which is consistent with ground state AFQMC for an open-shell system. Similarly, any degenerate partners of the conduction state are not included. After the Gram-Schmidt step, only the occupied orbitals are retained and included in the propagation and measurement, until the next time for orthogonalization when the ’s are re-inserted and the procedure repeated. The above procedure generalizes straightforwardly if the targeted excited state corresponds to “double” excitations or beyond.
We take simple trial wave functions directly from DFT calculations for the phaseless constraint Zhang and Krakauer (2003) and for the orthogonalization. The starting point for constructing at the selected point in the Brillouin zone (BZ) is the single Slater determinant DFT wave function for the corresponding ground state, . The orbitals are given by the eigenfunctions at based on a well-converged density integrated over -points. For a given excitation of spin- from an occupied orbital to an unoccupied orbital , we replace by . For insulators, the ground states are closed shell configurations, so this type of excited state Slater determinant will be spin-contaminated in general. To avoid this, we form a two-determinant singlet wave function:
where and are creation and destruction operators for unoccupied (conduction) and occupied (valence) states, respectively, and we assume that the spatial part of the valence and conduction orbitals are spin independent in . In degenerate cases, the trial wave function is constructed by considering all possible promotions among these orbitals, and the coefficients of each determinant are set equal. All our calculations are performed in the primitive cell. In supercells, the folding of bands creates additional mixing of crystal momentum eigenstates, whose decoupling will need further study.
Figure 1 illustrates our method in fcc silicon. The excited state corresponds to the excitation from band to the lowest lying conduction band at . The trial wave function was obtained from DFT with the local-density approximation (LDA), while for excited state is the singlet trial wave function of Eq. (2). For comparison, AFQMC results are also shown from free-projection, which does not impose the phaseless constraint Zhang and Krakauer (2003) and is exact for the ground state, but is eventually overwhelmed by the fermion sign problem. (Large runs with walkers were used in the free-projection calculations.) The ground state phaseless AFQMC results are in excellent agreement with the exact calculation, indicating very small systematic error from the phaseless approximation, consistent with earlier results Purwanto et al. (2009b). In contrast with the ground state, exactness is not ensured in free projection for the excited state; a more severe sign problem and onset of collapse to the ground state is seen. The phaseless and the additional orthogonality constraints stabilize the excited state calculation and yield an accurate excitation energy, which is within eV of the GW result Rohlfing et al. (1993) after correction of finite-size effects, as we discuss next.
Excitation energies from valence state to conduction state at any selected BZ point are calculated as the difference between the AFQMC total energy and the corresponding ground-state energy
Because the many-body calculations are performed in finite-size (FS) simulation cells, the energies have FS errors which must be removed Kent et al. (1999); Chiesa et al. (2006); Kwee et al. (2008); Ma et al. (2011); Drummond et al. (2008). FS corrections to the thermodynamic limit of an infinite-sized supercell can be obtained as post-processing corrections Kwee et al. (2008); Ma et al. (2011) from lower-levels of theory for both one-body (1b) and two-body (2b) effects. In the present excited state calculations, however, the creation of the electron-hole (eh) pair results in an additional bias, from the interaction between the particle and hole. We obtain the combined 1b and eh FS correction from DFT calculations:
where is the difference in band energies from a standard DFT calculation, using a dense -point grid, while is from self-consistent DFT calculations at paralleling the QMC calculations in Eq. (3). We correct the 2b FS error, which along with the 1b effect are substantially reduced here because is an energy difference between two states within the same simulation cell, using an LDA functional specially parametrized Kwee et al. (2008); Ma et al. (2011) for FS calculations
where the two excitation energies on the right are from LDA calculations paralleling the QMC, using the standard and the FS functionals, respectively. The sum of Eqs. (4) and (5) gives the total FS correction . The largest contribution is from , typically eV at most points in Si and diamond. Its largest value in Si is 0.35 eV at the point. In diamond, which has a large band gap, its largest value is 0.83 eV. The 2b correction is typically smaller; its largest value is 0.12 eV in Si and diamond, and 0.08 eV for the fundamental gap in ZnO.
We then obtain a quasiparticle band structure from a least-squares fit (Williamson et al., 1998) to the calculated many-body excitation energies :
where and run over the occupied (v) and unoccupied (c) states, respectively. The highest occupied quasiparticle energy is set equal to the corresponding DFT eigenvalue .
The full many-body Si band structure, with FS correction included, is shown in Fig. 2 and compared to LDA and , as well as to available DMC and experimental results Rohlfing et al. (1993); Williamson et al. (1998). Calculations for fcc Si were done at the experimental lattice constant of 5.431 . Both AFQMC and DFT calculations used a norm-conserving Kleinman-Bylander type separable nonlocal LDA pseudopotentials generated by the OPIUM code opi (), with a planewave cutoff Ry. Using a from LDA, AFQMC results correct the band gap problem and are in good agreement with experiment and with . The lowest band, which corresponds to the highest excitation energies, tends to be eV too low. For an imaginary-time projection method, its quality can be expected to decline for higher excited states. Also the simple singlet or the orthogonalization constraint may not be sufficient, as there are many states with similar energies.
The many body band structure of diamond, with FS correction included, is given in Fig. 3 together with available DMC and experimental results Towler et al. (2000); Rohlfing et al. (1993). The lattice constant of 3.567 was used. The calculations are similar to those in Si, except for a higher planewave cutoff Ry and the use of GGA pseudopotential and trial wave function. We have verified that the calculations are insensitive to the difference in at the level of LDA vs. GGA or a hybrid functional. The AFQMC results are again generally in very good agreement with and experiment.
Accurately calculating the fundamental band gap of wurtzite structure ZnO is very challenging. DFT LDA and GGA underestimate the gap by almost eV. Even the generally more accurate method underestimates the gap by - eV. There has been considerable discussion about the importance of various convergence issues, choice of pseudopotentials, and additional approximations such as the plasmon-pole model and the inclusion of self-consistency in the calculation Shih et al. (2010); Stankovski et al. (2011); Partoens et al. (2010). Results are also sensitive to the choice of pseudopotential, since there is large spatial overlap among , and wave functions of atomic Zn Partoens et al. (2010). To properly treat this, our Zn GGA pseudopotential was constructed with a Ne-core, thereby fully correlating the semicore , along with the states in the many-body calculation. Very conservative radial cutoffs of , , and bohr were used for , , and channels, respectively, resulting in a large planewave cutoff energy of Ry. The pseudopotential gives GGA optimized lattice parameters of and , and bulk modulus of GPa, all in good agreement with published GGA results Molepo and Joubert (2011).
Our AFQMC calculations were done at the above GGA optimized geometry. While hybrid functionals seem to perform better in ZnO, we chose to use simple trial wave functions from GGA in our AFQMC calculations to avoid any parameter tuning. The singlet form in Eq. (2) was used for the excited state. The main calculations were done with a time-step of . The Trotter error was then corrected by extrapolation from separate runs using a single-determinant . The calculated raw band gap was 2.54(14) eV. Since DFT severely underestimates the band gap, the FS correction is less straightforward in ZnO. At , for example, the single -point self-consistent GGA calculation would yield a metallic ground state. We extrapolated for within of the point, along two high symmetry lines. In doing so, we assume that the electron-hole effect does not introduce a discontinuity in the band dispersion, which is reasonable since it mainly relates to the simulation cell size. This yielded eV. (Had we used the LDA, which has a smaller gap, would be about 0.41(04) eV.) Ths 2b correction is well-behaved: eV. Adding the FS corrections yields our calculated band gap of eV.
The experimental equilibrium geometry is at and Desgreniers (1998). To compare our band gap (for the GGA equilibrium geometry) to experiment, we apply a correction of eV, which is the excitation energy difference given by GGA for the experimental and GGA-optimized geometries. (Hybrid B3LYP calculations give a correction of eV.) Table 1 compares our result with experimental values (3.30 Srikant and Clarke (1998), 3.44 Mang et al. (1995), and 3.57 eV Tsoi et al. (2006); Alawadhi et al. (2007)) and those from recent calculations Stankovski et al. (2011); Friedrich et al. (2011a); Berger et al. (2012); Betzinger et al. (2010); Uddin and Scuseria (2006). We note that, even with the small-core pseudopotential, there can be pseudoization and core-relaxation errors Gómez-Abal et al. (2008). The precise effect cannot be determined without further many-body calculations. However, recent studies Friedrich et al. (2011a); Berger et al. (2012) comparing all-electron and pseudopotential calculations indicate that the effect is approximately eV for a pseudopotential of similar quality to the one we have adopted. This would increase the AFQMC result for the fundamental band gap in ZnO to eV, within the range of experimental measurements Srikant and Clarke (1998); Mang et al. (1995); Tsoi et al. (2006); Alawadhi et al. (2007) listed in Table 1.
Various further improvements can be explored. We have used a planewave basis and norm-conserving pseudopotentials in these calculations. Other single-particle basis sets are straightforward to use, for example, with localized or natural orbitals. One could also work with a truncated set of orbitals from a lower-level of theory, or use a down-folded Hamiltonian directly to improve computational efficiency. The constraining virtual orbitals have been fixed in our calculations, but could potentially be allowed to dynamically evolve in some way in the propagation. It is reassuring that DFT trial wave functions have worked well. Clearly more elaborate multi-determinant wave functions can be used, for example, from diagonalization in a subspace formed by single and double excitations to conduction orbitals.
In summary, we have presented an AFQMC approach for the calculation of electronic excitations in solids, introducing an orbitally-based orthogonalization constraint in the phaseless AFQMC framework to stabilize the projection of excited states. Simple trial wave functions directly from DFT calculations were used for the constraint. Detailed many-body quasiparticle band structures can be calculated. In prototypical semiconductors (Si and diamond), the calculated band structures are in good agreement with those from calculations and with experiment. In the more strongly correlated and challenging wurtzite ZnO crystal, the calculated fundamental gap is in excellent agreement with the latest experimental data. The method is non-perturbative and free of empirical parameters, offering a possible path for general computations in correlated materials.
This work is supported by DOE (DE-FG02-09ER16046), NSF (DMR-1006217), and ONR (N000140811235; N000141211042). An award of computer time was provided by the Innovative and Novel Computational Impact on Theory and Experiment (INCITE) program, using resources of the Oak Ridge Leadership Computing Facility at the Oak Ridge National Laboratory, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC05-00OR22725. We also acknowledge the computing support from the Center for Piezoelectrics by Design. We thank W. Purwanto for many useful discussions, and Eric J. Walter for help with pseudopotentials.
- Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
- Kohn and Sham (1965) W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
- Perdew and Levy (1983) J. P. Perdew and M. Levy, Physical Review Letters 51, 1884 (1983).
- Becke (1993) A. D. Becke, J. Chem. Phys. 98, 1372 (1993).
- Onida and Rubio (2002) G. Onida and A. Rubio, Rev. Mod. Phys. 74, 601 (2002).
- Hedin (1965) L. Hedin, Phys. Rev. 139, A796 (1965).
- Shih et al. (2010) B.-C. Shih, Y. Xue, P. Zhang, M. L. Cohen, and S. G. Louie, Phys. Rev. Lett. 105, 146401 (2010).
- Stankovski et al. (2011) M. Stankovski, G. Antonius, D. Waroquiers, A. Miglio, H. Dixit, K. Sankaran, M. Giantomassi, X. Gonze, M. Côté, and G.-M. Rignanese, Phys. Rev. B 84, 241201 (2011).
- Partoens et al. (2010) H. D. Partoens, R. Saniz, D. Lamoen, and B, J. Phys.: Condens. Matter 22, 125505 (2010).
- Foulkes et al. (2001) W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal, Rev. Mod. Phys. 73, 33 (2001).
- Zhang and Krakauer (2003) S. Zhang and H. Krakauer, Phys. Rev. Lett. 90, 136401 (2003).
- Booth et al. (2009) G. H. Booth, A. J. W. Thom, and A. Alavi, The Journal of Chemical Physics 131, 54106 (2009).
- Schautz et al. (2004) F. Schautz, F. Buda, and C. Filippi, The Journal of Chemical Physics 121, 5836 (2004).
- Purwanto et al. (2009a) W. Purwanto, S. Zhang, and H. Krakauer, J. Chem. Phys. 130, 94107 (2009a).
- Williamson et al. (1998) A. J. Williamson, R. Q. Hood, R. J. Needs, and G. Rajagopal, Phys. Rev. B 57, 12140 (1998).
- Towler et al. (2000) M. D. Towler, R. Q. Hood, and R. J. Needs, Phys. Rev. B 62, 2330 (2000).
- Chang and Zhang (2008) C.-C. Chang and S. Zhang, Physical Review B 78, 165101 (2008).
- Suewattana et al. (2007) M. Suewattana, W. Purwanto, S. Zhang, H. Krakauer, and E. J. Walter, Phys. Rev. B 75, 245123 (2007).
- Kwee et al. (2008) H. Kwee, S. Zhang, and H. Krakauer, Phys. Rev. Lett. 100, 126404 (2008).
- Al-Saidi et al. (2006) W. A. Al-Saidi, H. Krakauer, and S. Zhang, Phys. Rev. B 73, 075103 (2006).
- Rohlfing et al. (1993) M. Rohlfing, P. Krüger, and J. Pollmann, Phys. Rev. B 48, 17791 (1993).
- Purwanto et al. (2009b) W. Purwanto, H. Krakauer, and S. Zhang, Physical Review B 80, 214116 (2009b).
- Kent et al. (1999) P. R. C. Kent, R. Q. Hood, A. J. Williamson, R. J. Needs, W. M. C. Foulkes, and G. Rajagopal, Physical Review B 59, 1917 (1999).
- Chiesa et al. (2006) S. Chiesa, D. M. Ceperley, R. M. Martin, and M. Holzmann, Physical Review Letters 97, 076404 (2006).
- Ma et al. (2011) F. Ma, S. Zhang, and H. Krakauer, Phys. Rev. B 84, 155130 (2011).
- Drummond et al. (2008) N. D. Drummond, R. J. Needs, A. Sorouri, and W. M. C. Foulkes, Physical Review B 78, 125106 (2008).
- (27) The OPIUM project, available at http://opium.sourceforge.net.
- Faleev et al. (2006) S. V. Faleev, M. van Schilfgaarde, T. Kotani, F. Léonard, and M. P. Desjarlais, Phys. Rev. B 74, 033101 (2006).
- Molepo and Joubert (2011) M. P. Molepo and D. P. Joubert, Phys. Rev. B 84, 094110 (2011).
- Desgreniers (1998) S. Desgreniers, Physical Review B 58, 14102 (1998).
- Srikant and Clarke (1998) V. Srikant and D. R. Clarke, Journal of Applied Physics 83, 5447 (1998).
- Mang et al. (1995) A. Mang, K. Reimann, and S. Rübenacke, Solid State Communications 94, 251 (1995).
- Tsoi et al. (2006) S. Tsoi, X. Lu, A. K. Ramdas, H. Alawadhi, M. Grimsditch, M. Cardona, and R. Lauck, Physical Review B 74, 165203 (2006).
- Alawadhi et al. (2007) H. Alawadhi, S. Tsoi, X. Lu, A. K. Ramdas, M. Grimsditch, M. Cardona, and R. Lauck, Physical Review B 75, 205207 (2007).
- Friedrich et al. (2011a) C. Friedrich, M. C. Müller, and S. Blügel, Physical Review B 83, 081101 (2011b); 84, 039906(E) (2011a).
- Berger et al. (2012) J. A. Berger, L. Reining, and F. Sottile, Physical Review B 85, 085126 (2012).
- Betzinger et al. (2010) M. Betzinger, C. Friedrich, and S. Blügel, Phys. Rev. B 81, 195117 (2010).
- Uddin and Scuseria (2006) J. Uddin and G. E. Scuseria, Phys. Rev. B 74, 245115 (2006).
- Gómez-Abal et al. (2008) R. Gómez-Abal, X. Li, M. Scheffler, and C. Ambrosch-Draxl, Physical Review Letters 101, 106404 (2008).