First we note that the rate of convergence is different for different frequencies even though a possible pattern has not been identified. This confirms the importance of monitoring the convergence in a frequency-dependent manner in order to avoid unnecessary iterations for those frequencies already below a given convergence threshold. Moreover the convergence of the simple fixed-point algorithm suggests that the rs-SCS problem is well-conditioned even for the dense cluster of atoms taken as test system. It is also of interest to probe the effects of the DIIS extrapolation coupled to the fixed-point iterative solution of the screening equations and this is shown in Figure . As expected, the DIIS extrapolation accelerates the convergence, in particular for the slow converging frequencies observed previously in Figure ??. The total number of iterations is reduced, for the full set of frequencies, by 55%, underlining the benefits of the DIIS inclusion. As mentioned in the theoretical section, an additional aspect considered to improve the performance of the algorithm is the re-use of purely geometric-dependent terms (referred to as G i j in Algorithm 1) in the on-the-fly building of T SR i j (iν n ): this allows for further decreasing the computational cost by up to 40% for a 10-points quadrature as shown in Figure . A final aspect concerning the efficient solution of the rs-SCS sets of equations is the choice of the proper cutoff radius from which the neighbor lists can be built. The short-range nature of T SR responsible for its sparsity can be seen in Figure where the isotropized rs-SCS polarizability of one of the argon atoms in a dimer is plotted as function of the inter-atomic distance. At inter-atomic distances above 3 Ångstroms the effects of the self-consistent screening becomes negligible suggesting that a cutoff radius of 4 Ångstroms is more than enough to include all the anisotropic effects arising from the neighboring atoms. We wish add to the discussion of the proposed iterative rs-SCS algorithm described in Algorithm 1 the comparison between its computational scaling and the one relative to the direct inversion approach in Eq.( ), Figure . The remarkable computational gain introduced by the iterative solution of the reduced rs-SCS Dyson-like equations shown in Figure is, however, only one of the advantages introduced: a further aspect to take into account is in fact the memory efficiency. The solution of Eq.( ) requires to store the upper-triangular part of Ā(iν) for all the set of frequencies employed, thus involving quadratic-scaling memory requirements with the system size. The variable reduction involved by Eq.( ), on the other hand, only requires 9N elements to be stored for each of the frequencies employed, making the needed memory allocation linear scaling with the system size as shown in Figure . The advantages introduced by the here proposed strategy thus concerns both the computational performance and memory allocation requirements.