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 GW and experimental results. For the challenging ZnO wurtzite structure, we obtain a fundamental band gap of 3.26(16) eV, consistent with experiments.

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 [17] and real material systems [11,[18][19][20].
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 3dbands interacting with 3s and 3p semicore states.
In AFQMC, we target an excited state of the manybody HamiltonianĤ with the imaginary-time projection |Ψ * (β) ∝ e −βĤ |ψ * T , where |ψ * T is a trial wave function of the excited state, and β ≡ N ∆τ with N the simulation time step and ∆τ the Trotter step size. The wave function at any imaginary-time can be thought of as |Ψ * (β) = φ c φ |φ ↑ ⊗ |φ ↓ . Each Slater determinant in the sum is a random walker, and their distribution gives a statistical representation of the coefficients c φ . Explicitly, |φ ↑ has the form (φ 1 , · · · , φ i , · · · , φ m ), where each φ i is an orbital (expressed in terms of the chosen single-partle basis) that evolves with β, and m 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 E * can be calculated by the mixed estimator [14]: which converges to the exact result at β → ∞. If |ψ * T 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) [14]. 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 i th valence orbital by a conduction orbital labeled by j (thus j > m). Each random walker is still an m-electron Slater determinant.
For the purpose of orthogonalization, however, we will regard it as the extended ordered list (φ 1 , · · · , φ i−1 ,φ i , φ i+1 , · · · , φ m ,φ m+1 , · · · ,φ j−1 , φ j ). Any orbital denoted by '-' is a virtual orbital whose only role is in the orthogonalization of the m occupied orbitals. As an approximation, we will use the corresponding orbitals in |ψ T , i.e.,φ i = φ T i . The choice of |ψ T is further discussed below. A Gram-Schmidt orthogonalization on this extended set ensures that the following orthogonality conditions are obeyed: i) φ i |φ j = 0; ii) φ i |φ k = 0 for k ∈ [i + 1, m]; iii) φ k |φ j = 0 for k ∈ [m + 1, j − 1]. If the valence stateφ i 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 φ j are not included. After the Gram-Schmidt step, only the m 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. . Phaseless AFQMC result for the ground state (red diamonds) is compared to exact free-projection (black circles). For the excited state, freeprojection (blue squares) collapses. Result from the new method with phaseless and virtual orthogonalization constraint is shown in red triangles, and compared to the GW excitation, indicated by violet dashed lines [21].
We take simple trial wave functions directly from DFT calculations for the phaseless constraint [11] and for the orthogonalization. The starting point for constructing |ψ * T (k) at the selected k point in the Brillouin zone (BZ) is the single Slater determinant DFT wave function for the corresponding ground state, |ψ T (k) . The orbitals φ T σ,i are given by the eigenfunctions at k based on a wellconverged density integrated over k-points. For a given excitation of spin-σ from an occupied orbital i to an unoccupied orbital j, we replace φ T i,σ by φ T j,σ . 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 a † j and a i 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 |ψ T . 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 i = 2 to the lowest lying conduction band j = 5 at k = (0.3, 0, 0). The trial wave function |ψ T was obtained from DFT with the local-density approximation (LDA), while for excited state |ψ * T 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 [11] and is exact for the ground state, but is eventually overwhelmed by the fermion sign problem. (Large runs with ∼ 2 × 10 6 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 [22]. 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 ∼ 0.3 eV of the GW result [21] after correction of finite-size effects, as we discuss next.
Excitation energies from valence state i to conduction state j at any selected BZ k 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 [19,[23][24][25][26]. FS corrections to the thermodynamic limit of an infinite-sized supercell can be obtained as post-processing corrections [19,25] 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 e g (k) ≡ e j (k) − e i (k) is the difference in band energies from a standard DFT calculation, using a dense k-point grid, while ∆E DFT ij (k) is from self-consistent DFT calculations at k paralleling the QMC calculations in Eq. (3). We correct the 2b FS error, which along with the 1b effect are substantially reduced here because ∆E QMC ij (k) is an energy difference between two states within the same simulation cell, using an LDA functional specially parametrized [19,25] 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 δE ij (k). The largest contribution is from δE eh,1b ij (k), typically ∼ 0.10 eV at most k 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 ǫ n (k) from a least-squares fit [15] to the calculated many-body where i and j run over the occupied (v) and unoccupied (c) states, respectively. The highest occupied quasiparticle energy is set equal to the corresponding DFT eigenvalue ǫ m (k) = e m (k). The full many-body Si band structure, with FS correction included, is shown in Fig. 2 and compared to LDA and GW , as well as to available DMC and experimental results [15,21]. 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 [27], with a planewave cutoff E cut = 25 Ry. Using a |ψ * T from LDA, AFQMC results correct the band gap problem and are in good agreement with experiment and with GW . The lowest band, which corresponds to the highest excitation energies, tends to be ∼ 1.5 eV too low. For an imaginary-time projection method, its quality can be expected to decline for higher excited states. Also the simple singlet |ψ * T 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 [16,21]. The lattice constant of 3.567Å was used. The calculations are similar to those in Si, except for a higher planewave cutoff E cut = 50 Ry and the use of GGA pseudopotential and trial wave function. We have verified that the calculations are insensitive to the difference in |ψ * T at the level of LDA vs. GGA or a hybrid functional. The AFQMC results are again generally in very good agreement with GW and experiment.
Accurately calculating the fundamental band gap of wurtzite structure ZnO is very challenging. DFT LDA and GGA underestimate the gap by almost 3 eV. Even the generally more accurate GW method underestimates the gap by 0.8 -1.1 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 GW calculation [7][8][9]. Results are also sensitive to the choice of pseudopotential, since there is large spatial overlap among 3s, 3p and 3d wave functions of atomic Zn [9]. To properly treat this, our Zn GGA pseudopotential was constructed with a Necore, thereby fully correlating the semicore 3s, 3p along with the 3d states in the many-body calculation. Very conservative radial cutoffs of 0.83, 1.02, and 1.13 bohr were used for s, p, and d channels, respectively, resulting in a large planewave cutoff energy of E cut = 180 Ry. The pseudopotential gives GGA optimized lattice parameters of a = 3.279 and c = 5.284Å, and bulk modulus of 128.7 GPa, all in good agreement with published GGA results [29].
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 |ψ * T in Eq. (2) was used for the excited state. The main calculations  [31][32][33][34]. Also shown are results from GGA, hybrid DFT [37,38], and GW [8,35,36]. were done with a time-step of ∆τ = 0.01Ry −1 . The Trotter error was then corrected by extrapolation from separate runs using a single-determinant |ψ * T . 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 k = Γ, for example, the single k-point self-consistent GGA calculation would yield a metallic ground state. We extrapolated δE eh,1b ij (k) for k within 0.1 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 δE 1b,eh The experimental equilibrium geometry is at a = 3.250 and c = 5.207Å [30]. To compare our band gap (for the GGA equilibrium geometry) to experiment, we apply a correction of +0.06 eV, which is the excitation energy difference given by GGA for the experimental and GGAoptimized geometries. (Hybrid B3LYP calculations give a correction of +0.10 eV.) Table I [32], and 3.57 eV [33,34]) and those from recent calculations [8,[35][36][37][38]. We note that, even with the small-core pseudopotential, there can be pseudoization and core-relaxation errors [39]. The precise effect cannot be determined without further many-body calculations. However, recent studies [35,36] comparing all-electron and pseudopotential GW calculations indicate that the effect is approximately +0.27 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 ∼ 3.53(16) eV, within the range of experimental measurements [31][32][33][34] listed in Table I.

GGA
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 func-tions have worked well. Clearly more elaborate multideterminant 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 GW 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.