Pairwise products of electron orbitals appear ubiquitously in quantum chemistry applications, such as the Fock exchange operator, the MP2 amplitude, and the polarizability operator. In a system with N_e electrons, the total number of orbital pairs is O(N_e^2) while the degree of freedom is only O(N_e). Many methods have been developed to compress the redundant information, among them the recently proposed interpolative separable density fitting (ISDF). In the context of hybrid density functional calculations, although ISDF significantly reduces the computation cost, it still requires an expensive QR factorization with column pivoting (QRCP) to select a set of interpolation points. In this work, we propose a new approach of finding interpolation points using centroidal Voronoi tessellation (CVT), which achieves comparable accuracy at much lower cost. Our method also enhances the smoothness of the potential energy surface in the context of ab initio molecular dynamics (AIMD) simulations with hybrid functionals.