A robust density fitting method for calculating Coulomb matrix elements over Bloch functions based on the 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 Hartree-Fock (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 (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 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.