A resolution-of-identity-based real-space implementation of the exact exchange algorithm was optimized to allow for much improved exploitation of sparsity and load balancing across ten thousands of parallel computational tasks. Results show drastically improved memory and runtime performance, scalability, and workload distribution on CPU clusters. The improvements pushed the simulation limits beyond 10,000 atoms, compared to an earlier implementation that reached system sizes around 1,000 heavy atoms . This new code implementation can perform computation of energy, forces, and stress for periodic and non-periodic systems for several fashions of hybrid density functionals. In addition, for materials including heavy elements, perturbative spin-orbit coupling can be combined with the hybrids . Due to the inherent O(N 3 ) scaling, the solution of the eigenvalue problem beyond 10,000 atoms becomes the bottleneck during the simulations.