Rather than using a grid-based algorithm to solve the TDKS equations, we follow the path established by Saalfrank, Luppi, and others, using atomcentered Gaussian basis sets to represent the density. This makes TDKS simulations accessible with widelyavailable, general purpose quantum chemistry codes, and the compactness of the Gaussian representation offers the possibility to extend TDKS calculations to larger molecules than are feasible in grid-based representations. Even for H 2 , where both approaches are feasible, the Gaussian representation is found to be at least an order of magnitude more efficient, and only that approach is scalable. This efficiency might be used, for example, to simulate HHG in a liquid or in a bulk crystal, both of which have been realized experimentally. For HHG in liquids, the role of molecular disorder has been emphasized, suggesting that model systems larger than a single molecule are required to understand this phenomenon.