To integrate over imaginary frequencies, we implemented a Gauss-Legendre (GL) quadrature scheme using 16 points, which has been previously shown to be effective at predicting accurate RPA energies. In this scheme, abscissa points in the standard GL interval x ∈ [-1, 1] are first mapped to the interval [0, 1] by the function y i = 1/2 + x i /2, where i is an index that runs over the number of imaginary frequencies used in the quadrature. This interval is then mapped to the desired integration range ω i ∈ [0, ∞[ by the function ω i = (-ln(1 -y i )) B , where B determines the density of frequency points at low energies and is either 2.5 for metals or 2.0 for other materials. The mapped ω i are normalized by (-ln(1 -max(y i )) B /E max , where E max determines the highest frequency point, and is set to 800 eV here. Due to the change in integration bounds, the quadrature weights are also transformed by the derivative of overall function mapping x i to ω i .