The overall computational cost of Eq.( ) arises from the AT SR p[k] σ matrix-vector products involved in the iterative procedure, where the diagonality of A allows the computation of AT SR on the fly at nearly no cost. The matrix-vector products described above can be performed very efficiently by exploiting the short-range nature of T SR , i.e. its sparsity. This is achieved by building specific short-range neighbor lists based on a cut-off radius which choice will be discussed and motivated in the numerical results section. A further aspect to consider in improving the matrix-vector products' efficiency is the way the T SR elements are built. We remind that Eq.( ) should be solved as many times as the number of frequencies (N ν ) needed in the GL quadrature employed to solve Eq. (14). In this regard it is possible to compute the purely geometry-dependent parts of T SR once and this allows their re-use in the on-the-fly building of frequency-dependent elements. This, as it will be analyzed later, decreases considerably the computational cost. Efficient matrix-vector products would be of little help if the overall convergence required a large number of iterations: it is therefore important to make use of a robust and efficient iterative scheme. The simple self-consistent solution of Eq.( ) can be considered as a fixed-point iterative procedure where the previous iteration's result is used as input for the next one. The rate of convergence of fixed-point iterative algorithms can be generally improved by constructing, at each iteration, a guess vector being a linear combination of previous iterations' resulting vectors where the coefficients of the linear combination are determined to minimize in a least squares sense the error. This extrapolation scheme known as diret inversion of the iterative subspace (DIIS) was introduced with the aim of accelerating and stabilizing the convergence of the Roothaan-Hall self-consistent field equations. The application of the DIIS together with fixed-point iterative methods in solving linear systems of equations makes the overall procedure equivalent to the Generalized Minimal Residual method (GMRES) which, belonging to the Krylov subspace class of iterative methods, is known for its fast and robust convergence properties. We mention that the coupled fixed-point/DIIS iterative schemes have proven successful in the context of modelling many-body polarization interactions and this motivates the use of this mixed procedure in the solution of Eq. (25). Finally, when the px , py and pz arrays are converged for all the N ν frequencies, the frequency-dependent screened atomic polarizabilities are isotropized (second equality in Eq.( )) and the set of { ωi } constituting the zero-point energy E 0 are obtained via Eqs.( ) and (2) respectively.