  We propose an efficient algorithm for density fitting of Bloch waves for
Hamiltonian operators with periodic potential. The algorithm is based on column
selection and random Fourier projection of the orbital functions. The
computational cost of the algorithm scales as $\mathcal{O}\bigl(N_{\text{grid}}
N^2 + N_{\text{grid}} NK \log (NK)\bigr)$, where $N_{\text{grid}}$ is number of
spatial grid points, $K$ is the number of sampling $k$-points in first
Brillouin zone, and $N$ is the number of bands under consideration. We validate
the algorithm by numerical examples in both two and three dimensions.
