  We generalise the non-affine theory of viscoelasticity for use with large,
well-sampled systems of arbitrary chemical complexity. Having in mind
predictions of mechanical and vibrational properties of amorphous systems with
atomistic resolution, we propose an extension of the Kernel Polynomial Method
(KPM) for the computation of the vibrational density of states (VDOS) and the
eigenmodes, including the $\Gamma$-correlator of the affine force-field, which
is a key ingredient of lattice-dynamic calculations of viscoelasticity. We show
that the results converge well to the solution obtained by direct
diagonalization (DD) of the Hessian (dynamical) matrix. As is well known, the
DD approach has prohibitively high computational requirements for systems with
$N=10^4$ atoms or larger. Instead, the KPM approach developed here allows one
to scale up lattice dynamic calculations of real materials up to $10^6$ atoms,
with a hugely more favorable (linear) scaling of computation time and memory
consumption with $N$.
