The Random Phase Approximation (RPA), which represents the fifth rung of accuracy in Density Functional Theory (DFT), is made practical for large systems. Energies of condensed phase systems containing thousands of explicitly correlated electrons and 1500 atoms can now be computed in minutes and less than 1 h, respectively. GPU acceleration is employed for dense and sparse linear algebra, while communication is minimized by a judicious data layout. The performance of the algorithms, implemented in the widely used CP2K simulation package, has been investigated on hybrid Cray XC30 and XK7 architectures, up to 16,384 nodes. Our results emphasize the importance of good network performance, in addition to the availability of GPUs and generous on node memory. A new level of predictivity has thus become available for routine application in Monte Carlo and molecular dynamics simulations.