Hybrid Density Functional Theory (DFT) has recently gained popularity as an accurate model of electronic interactions in chemistry and materials science applications. The most computationally expensive part of hybrid DFT simulations is the calculation of exchange integrals between pairs of electrons. We present strategies to achieve improved load balancing and scalability for the parallel computation of these integrals. First, we develop a cost model for the calculation, and utilize random search algorithms to optimize the data distribution and calculation schedule. Second, we further improve performance using partial data-replication to increase data availability across cores. We demonstrate these improvements using an implementation in the Qbox Density Functional Theory code on the Mira Blue Gene/Q computer at Argonne National Laboratory. We perform calculations in the range of 8k to 128k cores on two representative simulation samples from materials science and chemistry applications: liquid water and a metal-water interface.