Density fitting in periodic systems: application to TDHF in diamond and oxides

A robust density fitting method for calculating Coulomb matrix elements over Bloch functions based on calculation of two- and three-center matrix elements of the Ewald potential is described and implemented in a Gaussian orbital basis in the Exciton code. The method is tested by comparing Coulomb and exchange energies from density fitting to corresponding energies from SCF HF calculations for diamond, magnesium oxide and bulk Ne. Density fitting coefficients from the robust method are compared to coefficients from a variational method applied to wave function orbital products in bulk Ne. Four center Coulomb matrix elements from density fitting are applied to time dependent Hartree-Fock (TDHF) calculations in diamond, magnesium oxide and anatase and rutile polytypes of titanium dioxide. Shifting virtual states downwards uniformly relative to occupied states and scaling the electron-hole attraction term in the TDHF Hamiltonian by 0.4 yields good agreement with either experiment and/or Bethe-Salpeter equation calculations. This approach mirrors similar 'scissors' adjustments of occupied and virtual states and introduction of a scaled electron-hole attraction term in some time dependent DFT calculations.


I. INTRODUCTION
Excitons in condensed matter arise from a balance of electron-hole Coulombic attraction and electron-hole pair hopping. The former is mediated by finite wavevector q ladder matrix elements over four Bloch functions and, for q → 0 excitons which couple to light, the latter is mediated by ring matrix elements. Computation of these matrix elements is expensive and thus far has mainly been done in condensed phases using plane wave basis sets. Here we describe a method for computing these matrix elements using a Gaussian orbital basis and apply them to time dependent Hartree-Fock (TDHF) calculations in diamond, magnesium oxide and titanium dioxide. When a local orbital basis set is used for a periodic system, methods capable of treating the long range nature of the Coulomb interaction between local orbital basis functions are essential. Furthermore, four center matrix elements over the local orbital basis must be transformed to a Bloch orbital basis, which may also be expensive. Density fitting (DF) of products of wave function orbitals is a long established technique in both finite 1-10 and periodic systems [11][12][13][14][15][16][17][18][19][20][21][22][23][24][25] which can account for the long range Coulomb interaction and reduces the time required for integral calculation and transformation.
DF factorizes calculation of four center matrix elements of the Coulomb potential into products of two center Coulomb and three centre Coulomb or overlap integrals, depending on metric choice. This reduces the transformation from the local orbital basis to the Bloch orbital basis from a four orbital to a two orbital problem. When this is combined with space group symmetries of pairs of k points (k and k + q) in reciprocal space, high density k point meshes (up to 14 x 14 x 14 in this work) can be reached. All calculations reported were done using the Exciton code 26,27 which uses a local Gaussian orbital basis for SCF Hartree a) Electronic mail: Charles.Patterson@tcd.ie Fock (HF) and time dependent Hartree Fock (TDHF) calculations on finite and periodic systems. The code is also capable of performing GW and Bethe-Salpeter equation (BSE) calculations in finite systems 27 .
Here we report periodic TDHF calculations on diamond and three oxide compounds. First of all the theory underlying the methods used is described, secondly the accuracy of the DF approach used is evaluated by comparing Coulomb and exchange energies derived from DF to the corresponding SCF energies. The method used is a robust DF method in which the error in the difference between the fitted and true densities of products of pairs of orbitals is minimized with a Coulomb metric. However, net charges associated with these products are not reproduced exactly. The method of Lagrange multipliers is used to constrain fitted densities to their exact values in a variation DF method. Differences in constrained and unconstrained fitting coefficients are evaluated for bulk f cc Ne. In the following section, results of TDHF calculations for the systems just mentioned are reported. It is well known that differences in energy eigenvalues from DFT or HF calculations significantly under or overestimate quasi-particle band gaps in solids. This can be corrected in a GW calculation. However, since no periodic GW method is yet available for the Exciton code, an approach which has been successfully employed in TDDFT calculations 28 is adopted. Virtual states are shifted down by a constant amount and electron-hole attraction matrix elements are scaled by a factor of 0.4 for both diamond and the oxides studied. This is similar to the approach of Botti et al. who introduced a term −α/q 2 into the TDDFT exchangecorrelation kernel for a range of insulating materials 28 . Scaling of these matrix elements replaces screening by the static, inverse dielectric function 29 .
The charge densities which are fitted in this work are products of crystal orbitals, which are Gaussian orbitals with Bloch translational symmetry. Our DF approach employs an auxiliary basis in which HF Bloch functions are expanded. An alternative approach is to use the properties of spherical harmonics to reduce the number of orbital products 30 . Since individual Gaussian orbitals have limited extent, the net charge associated with a particular orbital product is finite and is equal to a Fourier transformed overlap matrix element for that pair of Gaussian orbitals.
In general these charge distributions have nonzero charge and dipole moment. The importance of treating charge distributions in periodic systems with zero monopole, dipole and quadrupole moment 31,32 has frequently been emphasized in work on DF in periodic systems 24,25,33 . Here we adopt an approach to fitting these orbital products in which this is taken into consideration using a conventional Ewald approach 32 . Ring or bubble matrix elements in the TDHF Hamiltonian are constructed from orbital products at single k points and are fitted via unmodulated Ewald sums (q = 0). Ladder matrix elements (i.e. the electron-hole attraction) are constructed from orbital products at k points separated by wave vector q and are fitted via modulated Ewald sums. In the former the G = 0 term is omitted and agreement between DF and SCF Coulomb energies is of the order of 20-40 µH per atom. In the latter the q → 0 limit must be taken. When the TDHF Hamiltonian is set up for a periodic system, this limit can only be reached by extrapolation.
Differences in Coulomb and exchange energies from DF and from SCF energies in this work are similar to those reported previously. For example, Burow et al. reported DF calculations of Coulomb energies in molecules, molecular crystals, graphite and diamond 18 . They found differences ranging from 4 to 37 µH per atom in finite, molecular systems and from 3 to 51 µH per atom in periodic systems. Sun et al. 25 calculated the Coulomb and exchange energies of a H lattice with the cubic diamond structure. Using a Gaussian orbital auxiliary basis they reproduced the SCF Coulomb energy to within 100µH and the exchange energy to 1mH per atom. For Si in the diamond structure they reproduced the HF energy per cell to within 100 µH. Supplementing the auxiliary bases with plane waves allowed energies to be reproduced to much greater accuracy (nH for Si). Milko et al. 15 reported differences in HF energies etc. in trans-polyacetylene (t-PA) and similar 1-D polymers of 60 µH per unit cell and minimal errors in the HF band gap less than 1 meV.
Computation of Coulomb and exchange matrix elements using DF for post SCF methods such as Möller-Plesset methods has been reported by Katouda and Nagase 20 for trans-polyacetylene using Gaussian and Poisson DF basis sets. They found correlation energy differences associated with DF as small as 30 µH and ranging up to 3 mH using mixed Gaussian and Poisson auxiliary basis sets. Maschio results from minimization of the quadratic error, where v(r, r ′ ) is the Coulomb potential, ∆ρ(r) = ρ(r) −ρ(r) is the difference in true and fitted densities and w(r − r ′ ) = v(r, r ′ ) is the Coulomb metric. Electron-hole hopping and electron-hole attraction matrix elements needed for TDHF calculations are, and respectively. Ψ vk (r) and Ψ ck+q (r) are valence and conduction band states at wave vectors k and k+q. Bloch functions, Ψ vk (r), are expanded as linear combinations of phase modulated local orbitals, φ m (r − R)e ik.R , with lattice translation vector, R, which are referred to as crystal orbitals (CO) with wave vector, k, and expansion coefficients, d vk m . Expansion of Bloch function products in Eq. 3 and 4 leads to CO products, the charge densities which are fitted in this work. An auxiliary basis of the form, where the auxiliary basis, χ q α (r), is a more extensive Gaussian orbital basis set than the CO basis in Eq. 5 is used to expand the Bloch functions.
The condition that auxiliary function expansion coefficients, c k,q α , minimize the error, in the Coulomb integral ρ k,q mn (r)|ρ * k,q rs (r ′ ) is, Manipulation of the lattice sums in Eq. 6, 7 and Eq. 8 using lattice translational invariance turns Eq. 9 into the following lattice modulated Ewald sums 32,35 (Appendix A), where the periodic Ewald potential is shown explicitly for clarity. The potential, is a lattice-modulated generalization of the familiar Ewald potential 36 , where γ is the splitting parameter in Ewald's method and Ω is the unit cell volume. Defining, the equation which yields the coefficients which minimize the error in Eq. 8 becomes c k,q α V q αβ = V k,q mnβ and the expansion coefficients are, where V q−1 αβ is the matrix inverse of V q αβ . Substitution of the expansion of the density, ρ k,q mn (r) = c k,q α χ q α (r), yields

B. Variational Density Fitting
For some applications, such as use of DF in self-consistent field calculations, it is desirable or essential that fitted densities not only minimize errors in Coulomb and exchange energies, but also conserve charge. The robust fit method obtained by minimizing the error in Eq. 9 is not variational 2,4 in that while errors in the electrostatic energy in Eq. 1 are minimized, the fitted charge densities are not constrained to equal the densities being fitted. The constraint that the integrated charge densities Q k,q mn = drρ k,q mn (r) and Q k,q mn = c k,q α drχ q α (r) be equal is imposed by minimizing the functional, Where λ and µ are Lagrange multipliers. The constraint, reduces to, using lattice translational invariance (Appendix B). The term on the left in Eq. 19 is the Fourier transform of the overlap matrix at wave vector k + q, S k+q mn . The bar on the coefficient c k,q α indicates a constrained coefficient. The resulting linear equations, which replace Eq. 14, are, where α = drχ α (r) is nonzero for l = 0 auxiliary functions and zero otherwise. For an auxiliary basis function with a single, normalized Gaussian this equals (2π/a) 3/4 , where a is the Gaussian exponent. Eq. 20 must be solved for each unique k and k + q pair, however, decomposition of the matrix on the left before solution of these equations needs to be performed only once for each unique q vector. Here we investigate how well the unconstrained, robust density fitting of Eq. 14 satisfies charge conservation for bulk fcc Ne by comparing robust and variational DF expansion coefficients. All other calculations reported here used robust, non-variational fitting.

C. Exchange and Coulomb energies
Comparison of Coulomb and exchange energies from SCF and DF calculations presents a useful means of evaluating errors in density fitted electrostatic energies in periodic systems. Coulomb and exchange energies, E C and E x , are obtained by replacing conduction band states by valence states in Eq. 3 and 4 and summing over wave vectors, k and k + q, and In order to transform the density fitted matrix elements over the local orbital basis in Eq. 16 into matrix elements over Bloch functions, where d vk m and d v ′ k+q n are Bloch function expansion coefficients and the superscript q indicates the q vector used in V k,q αmn . Coulomb and exchange energies in Eq. 21 and 22 become, and Differences in total Coulomb and exchange energies obtained using these expressions and from SCF calculations used to generate Bloch functions in Eq. 24 and 25 are reported below. These are in the 25 -50µH per atom range for Coulomb energies and 1mH per atom range for exchange energies extrapolated to infinite q sampling density.

D. Small q limit
The divergent nature of the V q αβ and V k,q mnβ matrix elements in Eq. 12 and 13 around q = 0 requires special attention. The contribution to the Coulomb energy at q = 0 is straightforward. The first term on the right in Eq. 11 with G = 0 is replaced by −π/γΩ in a 3-D periodic system 36 . Choosing γ → ∞ makes the real space sum on the right tend to zero, leaving the Fourier expansion of the Coulomb potential on the right.
Expansion of the complex exponential in this term for small q and G = 0 as e iq.(r−r ′ ) ≈ 1+iq.(r−r ′ ) allows small q contributions for various interactions to be determined. Valence and conduction states in Ψ * vk and Ψ ck or Ψ * v ′ k+q and Ψ c ′ k+q in Eq. 3 (ring diagrams) are orthogonal at small q. Hence the leading contribution in this case is from iq.(r − r ′ ). This term (of order meV) leads to splitting of longitudinal and transverse excitons 29 and is omitted in this work. For valence states Ψ * vk and Ψ v ′ k+q which occur in the exchange energy (Eq. 22) and electron-hole attraction (ladder diagrams, Eq. 4) and conduction states Ψ * c ′ k+q and Ψ ck in ladder diagrams, identical states (v = v ′ and c = c ′ ) have unit overlap as q → 0 while distinct states are orthogonal. The former have a contribution from the leading term in the expansion of e iq.(r−r ′ ) and the latter in the iq.(r − r ′ ) term.
The q → 0 limit of the exchange energy and the electron-hole attraction term in the TDHF Hamiltonian (Eq. 4) is obtained by assuming that matrix elements around q = 0 vary slowly with q and may be approximated by the average value of the analytically integrated Coulomb potential around q = 0. Averaging 4π/Ωq 2 over a sphere with volume Ω BZ /N 3 yields, N 3 is the number of q points in a regular Monkhorst-Pack net 37 . This energy is used for each state where v = v ′ and q → 0 in Eq. 22. There are methods for calculating the total exchange energy more accurately than simply by sampling a divergent function of q. Gygi and Baldereschi 38 introduced a divergent, periodic function which is added and subtracted from the divergent small q limit of the Coulomb potential to produce a smoothly varying potential which is integrated numerically and a divergent, periodic term which is integrated analytically over the Brillouin zone. Sorouri et al. 39 , Carrier et al. 40 and Duchemin and Gygi 41 later adopted a simpler auxiliary function, which may be applied to all crystal systems. Spencer and Alavi 42 adopted a function which has a spherical real space cutoff.
Here, however, the method of calculating TDHF spectra requires sampling of the electronhole Hamiltonian on a regular Monkhorst-Pack net. The sampling method described above is used to calculate the contribution to the exchange energy and the electron-hole Hamiltonian. Calculation of the exchange energy for a series of nets of increasing k point density allows the exchange energy to be extrapolated to infinite density. This approach is used in results reported in Section V. Extrapolated exchange energies are compared to the exchange energy from the SCF calculation, which generated the Bloch functions used in density fitted integrals.
The exchange energy in these SCF calculations is obtained from four-center, real-space integrals and the real-space density matrix, P CA jl is the jl element of the real-space density matrix at lattice vector C−A. This approach relies on the convergence of the real-space density matrix with lattice vector range, |CA| and |B|. The real-space density matrix is exponentially localized for a gapped material and the localization length reduces with increase of the band gap. In this work only wide gap materials are studied and the SCF exchange energy is well converged for these systems.
It is worth noting that the Ewald method of calculating exchange energies (Eq. 22) does not rely on convergence of the density matrix in real space because the Ewald method sums interactions to infinite range. It may therefore be a superior method of calculating the HF exchange operator in SCF calculations on narrow gap or conducting systems. Indeed, products of Bloch functions in matrix elements in many-body calculations are delocalized over all space, and therefore calculation of converged electron-hole attraction matrix elements in crystalline systems is not possible using real-space four-center integrals of the kind in Eq.
28. An Ewald approach such as that used here is essential. Casting the Coulomb interaction into reciprocal space in Eq. 11 by choosing γ → ∞ has the disadvantage that a very large number of G vectors must be used in order to obtain convergence for even relatively low Gaussian exponents, (of order 1, say). Hence, the mixed real and reciprocal space method advocated here has a number of important advantages.

E. Time dependent Hartree-Fock method
The TDHF equations are usually expressed as the following generalized eigenvalue problem, Here the Tamm-Dancoff approximation (TDA) to the TDHF equations is used throughout in which the off-diagonal B blocks are omitted. In this case the TDHF problem reduces to a standard eigenvalue problem with eigenvectors X and eigenvalues Λ. The matrix elements which appear in the A block of the TDHF Hamiltonian matrix are given in Eq. 3 and 4. The A matrix also contains differences in single particle eigenvalues on the diagonal. Results are reported in Section V where either differences in HF eigenvalues are used or where these differences are reduced by a constant 'scissors' shift of around 10 eV.

III. COMPUTATIONAL METHODS
The implementation of this method in the Exciton code is described in this section. Firstly, the method of calculation of finite-q matrix elements which appear in ladder diagrams is given, followed by zero-q matrix elements which appear in ring diagrams. Diagonalization of the A matrix in Eq. 29 is performed using the pzheevx routine in Scalapack 43 . Scalapack uses a block-cyclic matrix distribution. Here the row and column blocksize in this distribution is chosen to be the number of transitions per k-point, i.e. the number of occupied times the number of virtual states per k point in the active space. This means that the matrix to be diagonalized is split over cores by q point and all matrix elements corresponding to a particular (k, k + q) pair are sent to the same core.
(1) HF-SCF calculations are performed using a set of unique k-points in an N x N x N Monkhorst-Pack net. Wave functions at symmetry equivalent k-points may be generated by rotation of these unique k-points using one known symmetry operator for each equivalent k-point. Phases of wave functions at unique and symmetry equivalent k-points are therefore unique for each wave function and k-point.
(2) Lists of unique (k, k + q) pairs are generated using the little group of q, whose symmetry operations leave q invariant. A set of symmetry-unique q points in the Brillouin zone is identified, then k points which are unique under the little group for that q vector are identified. Finally, all symmetry equivalent (k, k + q) pairs are generated using the full set of space group operations. Unique instances of these pairs are stored along with the symmetry operators which generated them from a given unique (k, k + q) pair.
(3) Matrix elements over atomic orbitals in Eq. 12 and 13 for unique (k, k + q) pairs are calculated. These are transformed to matrix elements over Bloch functions at these points by multiplying in the d vk m expansion coefficients (Eq. 5). (4) Computation of matrix elements V k,q mnβ is distributed over cores by unique q point. V k,q mnβ matrix elements for unique (k, k + q) pairs are computed and these are rotated to all symmetry equivalent equivalent (k ′ , k ′ + q ′ ) pairs. Each wave function in each (k, k + q) pair must be rotated twice; once to a unique (k, k + q) pair and once more from there to an equivalent (k ′ , k ′ + q ′ ) pair.
Two consequences of this are that (i) wave functions generated by products of symmetry operations may differ by a phase factor from the definitive phase in (1) and (ii) degenerate wave functions generated by rotation by the single, unique operator in (1) may be mixed when rotated twice. Both (i) and (ii) are accounted for by obtaining the overlap matrix of wave functions rotated by the single, unique operator with those rotated twice. Columns of this overlap matrix yield the linear combination of matrix elements which must be used to ensure the uniqueness of each symmetry equivalent wave function product at k ′ and k ′ + q ′ .
(5) Once a set of V k,q mnβ matrix elements and its symmetry equivalents have been computed, they are distributed to cores according to the Scalapack block-cyclic distribution mentioned above using an MPI-Send blocking send.
(6) Calculation of q = 0 matrix elements which occur in ring diagrams is simplified as wave functions at only one k point occur in each matrix element as q = 0. Matrix elements calculated at the unique set of k points mentioned in (1) can therefore be used to generate all equivalent (k ′ , k ′ ) pairs in ring diagrams.

IV. BASIS SETS
The wave function basis sets used in ths work are the DEF2-TZVP basis sets of Weigend and Ahlrichs 44 , modified for the solid state by removing diffuse basis functions or increasing their exponents and the auxiliary basis sets are the DEF2-TZVP-RIFIT basis sets 45 . Basis functions with angular momenta with l values greater than l = 4 (h functions and higher) were omitted. The DEF2-SVP, -TZVP and -QZVP basis sets for Ne used in Section V B were used without modification. Modifications to DEF2-TZVP basis sets for C, O, Mg and Ti are given in the Supplementary Information.

V. RESULTS
In this Section differences in Coulomb and exchange energies obtained from SCF calculations and by DF are compared for diamond, magnesium oxide and bulk fcc Ne. Fitted charges, Q k,q mn and Q k,q mn (Eq. 18), obtained by solving Eq. 14 and 20 are computed for bulk fcc Ne. The latter charges are equal to their exact numerical values, S k+q mn , owing to the constraints imposed in Eq. 20. Finally, results of TDHF calculations for these systems and two polymorphs of TiO 2 are presented and analyzed.

A. Fitted density Coulomb and exchange energies
Differences in Coulomb and exchange energies obtained from SCF calculations and by DF are shown in Fig. 1 as a function of Brillouin zone sampling frequency, 1/N . Differences in Coulomb energies are expected to be roughly independent of sampling frequency (Section II D) as they depend only on matrix elements at q = 0. Fig. 1 shows that differences in exchange energies from SCF calculations and DF scale with sampling frequency. Extrapolation of the DF exchange energy is needed for comparison to the SCF exchange energy.
SCF and DF energy differences in Coulomb energies obtained from SCF calculations and DF are 45, 55 and 56 µH, respectively, for diamond, MgO (around 20-25 µH per atom) and bulk Ne. These values are comparable to those obtained for molecular density fitting 6 using a series of attenuated Coulomb metrics. Extrapolated differences in SCF and DF exchange energies are -0.52 mHa for diamond, 3.26 mHa for MgO and 0.25 mHa for bulk Ne.

B. Charge conservation in fitted densities
The total charge associated with the density, ρ k,q mn (r), defined in Eq. 6 is the overlap matrix element at wave vector k + q, S k+q mn (Eq. 19). Eqs. 14 and 20 were solved for unconstrained coefficients, c k,q α , and constrained coefficients, c k,q α and Lagrange multipliers, λ k,q , for bulk Ne at (k, q) = (0, 0). Unconstrained coefficients are compared to overlap matrix elements, S 0 mn , in Fig. 2 for DEF2-SVP,-TZVP and -QZVP basis sets 45 and their corresponding auxiliary basis sets. Orbital combinations which produce a monopole on site are < l|l > and < l|l + 2n > where l is the orbital angular momentum and n is a positive integer (a product of two 2p x states or a 2p x state times 4f x(x 2 −3y 2 ) state on a given site has a net monopole, for example). < s|d >, l = 0, n = 1 differences are all below 10 −12 and are not shown. Absolute differences are shown in Fig. 2  from a SVP tp QZVP basis; although the auxiliary basis improves from SVP to QZVP, the number of wave function basis functions being fitted increases significantly.

C. Dielectric functions from TDHF
It is well known that screening of the exchange interaction in the GW approximation, for example, is necessary to predict the particle-hole gap correctly in solids and that the bare electron-hole attraction term in optical excitations in materials is screened in BSE calculations 29 . The A and B matrices in the TDHF method are equivalent to those in a GW /BSE calculation except that: particle-hole energy differences on the diagonal of the A matrix are HF eigenvalues in TDHF, rather than GW quasi-particle energies; the electronhole interaction is the bare interaction in TDHF whereas it is screened in BSE. In this work, screening of the electron-hole attraction term is sufficiently important, even in wide gap insulators such as diamond and MgO, that the experimental spectrum can only be recovered if it is reduced significantly. If it is not scaled in diamond or MgO the lowest excitation is a strongly localized Frenkel exciton, split off from the continuous part of the dielectric function rather than a Wannier exciton with absorption enhancement at lower frequencies in the dielectric function. The former typically occurs in systems such as rare gas solids, but not in bulk oxides or semiconductors.
Electron-hole attraction matrix elements were uniformly scaled by a factor of 0.4 (i.e. with no q dependence) for diamond, MgO and two polytypes of TiO 2 . This is similar to the approach of Botti et al., who introduced a term −α/q 2 into the TDDFT exchange-correlation kernel for a range of insulating materials 28 . In that work, for diamond the optimal value of α was 0.6 and for MgO it was 1.8. Optimal values for six tetrahedral semiconductors and MgO scaled roughly with 1/ǫ ∞ 28 . Here, a constant scaling by α = 0.4 was found to give satisfactory agreement with experiment for diamond and the oxides studied, provided a suitable shift of virtual states was used. HF band structures of the materials studied in this work are shown in Fig. 3 and scaling factors, α, shifts of virtual states and lowest TDHF excitation energies, with and without scaling and shifting, are summarized in Table  I. For diamond and MgO spectra were averaged over several calculations with high density  The TDHF spectrum of diamond obtained after shifting virtual states and scaling electron-hole matrix elements is shown in Fig. 4. Four valence and the lowest four conduction band states were used in the TDHF calculation. The overall width of spectral features is greater than experiment by around 2.2 eV. The lowest energy excitation in a TDHF calculation with no shift of virtual states or scaling of the electron-hole attraction matrix elements is 12.34 eV and with scaling described below it is 6.56 eV.
The position of the main peak in the dielectric function occurs at 13.9 eV in the TDHF calculation, while it occurs at 11.8 eV in experiment. The Γ point HF band gap is 14.93 eV and is combined with a downward shift in virtual states of 7.7 eV, so that the particle-hole gap before the TDHF calculation is 7.23 eV. This may be compared with the G 0 W 0 quasiparticle gap at Γ in diamond of 7.5 eV 47 , which used a DFT-LDA band structure as input.
In that case the width of the spectral features is underestimated and the main peak in the dielectric function is underestimated by around 1 eV. The trends of over and underestimating widths of spectral features is likely to be the result of similar over and underestimation of valence and conduction band widths by HF and DFT-LDA approximations. The HF valence band width of diamond has previously been reported to be 28.67 eV 48   The TDHF dielectric function of MgO after shifting virtual states and scaling matrix elements is shown in Fig. 5. Four valence and four conduction band states were used in the TDHF calculation. The Γ point HF band gap is 19.38 eV and when combined with a virtual state shift of 11.0 eV, results in a particle-hole gap of 8.38 eV before the TDHF calculation. This compares with a converged G 0 W 0 quasi-particle gap of 7.9 eV 51 . The experimental optical gap of MgO is 7.83 eV 50,52 and the lowest energy TDHF excitation occurs at 7.59 eV with the shifted particle-hole gap. If no virtual state shift or matrix element scaling is used, the lowest energy excitation is a strongly localized exciton at 14.71 eV ( Table I). The main features of the dielectric function are reproduced in the TDHF calculation, including strong enhancement of optical absorption at the threshold energy and three higher energy peaks in the range to 25 eV, which are also reproduced in a BSE calculation 28 .

Anatase and Rutile TiO2
The TDHF dielectric functions of anatase and rutile polytypes of TiO 2 are compared to experiment in Figs. 6 and 7. Calculations were performed using a 6 x 6 x 6 Monkhorst-Pack net. All 12 predominantly O 2p valence states and 10 predominantly Ti 3d conduction band states were included in TDHF calculations and the virtual state shift was 8.2 eV for both polytypes.
The Γ point HF band gap for anatase is 13.71 eV and the indirect gap at the X and Γ points is 13.25 eV. When combined with a virtual state shift of 8.2 eV, the direct particlehole gap at Γ is 5.51 eV. This may be compared to values of direct G 0 W 0 quasiparticle gaps at Γ for anatase TiO 2 of 3.45 56 , 4.14 57 and 4.29 eV 58 and indirect gaps of 3.73 56 , 3.56 57 and 3.83 eV 58 using plane wave basis sets and a HSE06 hybrid density functional 56 , DFT-LDA 57 or DFT-PBE 58 density functionals. For rutile the HF gap is 13.02 eV and the particle-hole gap after the virtual state shift is 4.82 eV. This may be compared to values of G 0 W 0 quasiparticle gaps of 3.46 eV 56 , 3.38 eV 57 and 3.59 eV 58 . There is excellent agreement between the TDHF calculations and experiment 53 for both polymorphs, in terms of peak position and strength. The overall shapes of the spectra with the incident electric vector parallel or perpendicular to the c axis are reproduced.
Two sets of experimental data have been included in Fig 7 for rutile TiO 2 54,55 . They show the low energy peak in the experimental dielectric function with the electric vector parallel to the c axis around 4 eV rising to 15 to 18. The TDHF calculation shows a somewhat narrower peak rising to over 30. As noted by Landmann et al. 56 , the earlier dielectric function measurement for rutile (Ref. [55], Expt II) shows only a shoulder between 5 and 10 eV for the electric field vector parallel to the c axis, while TDHF (this work), BSE 56 and a more recent experiment (Ref. [54], Expt I) 54 show a second peak in that range, although these experimental data do not extend above 8 eV. This second peak is also found in the BSE calculations mentioned above 56,57,59 . A DF method for calculating ring and ladder Coulomb matrix elements over extended Bloch functions has been presented. Charge densities which arise from products of crystal orbitals at k and k + q are expanded in an auxiliary basis of Gaussian crystal orbitals. Fitting of these densities requires calculation of two and three center matrix elements of these crystal orbitals over a lattice modulated Ewald potential. In a robust fit of the densities, expansion coefficients are obtained by inverting the two center matrix elements. Alternatively, a variational fit of these densities, in which the integrated charge densities are constrained to have their exact values, requires solution of a set of linear equations with Lagrange multipliers.
We have demonstrated that the DF method described is able to recover Coulomb and exchange energies for several periodic systems to a similar level of accuracy as has been reported elsewhere for periodic systems 15,18,25 . For light atoms such as C, O or Mg the SCF Coulomb energy is reproduced to within 50 µH per atom, including core contributions. The exchange energy extrapolated to infinite sampling density agrees with the SCF exchange energy to within 1 mH per atom.
TDHF calculations presented here use uniform scaling of the electron-hole attraction matrix elements rather than a q-dependent inverse dielectric matrix and shift the HF virtual states downward by a constant amount in order to achieve agreement with experiment for diamond and three oxide compounds. The shifted virtual state band gaps used are similar to G 0 W 0 band gaps, but may exceed the G 0 W 0 band gap. There may therefore be a range of scaling of the electron-hole attraction matrix elements and virtual state shift which yield comparable agreement with experiment. For MgO, a particle-hole gap around 0.4 eV greater than the G 0 W 0 gap and gave good agreement with experiment for the onset of optical absorption. For anatase and rutile TiO 2 , values around 1.5 eV greater than G 0 W 0 values were used and good agreement with onset of optical absorption was predicted.

VII. DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.
Making the substitutions, B = A + A ′ , C = A + C ′ and r ′ − C ′ = r ′′ and using lattice translational invariance, the left hand side becomes, as given in Eqn. 12 and 13.
Appendix B: Derivation of Eq. 19 The constraint on charge densities of orbital products in Eq. 18 may be expressed as, dr(ρ k,q mn (r) − c k,q α χ q α (r)) = 0, Substituting for ρ k,q mn (r) and χ q α (r) using Eq. 6 and 7 and replacing the lattice vector B by A + B ′ yields, or, where, and,