Accurate descriptions of intermolecular interactions are of great importance in simulations of molecular liquids. We present a novel method that combines in a variational formulation the accuracy of the Harris functional approach with the computational efficiency of linear-scaling subsystem density functional theory (DFT). While subsystem DFT methods allow for efficient large-scale simulations, they often lack in accuracy. The proposed method allows for simulations with millihartree accuracy in comparison to the Kohn-Sham DFT reference, at near linear-scaling efficiency. Embedded in the CP2K program package, the method is designed to enable ab initio molecular dynamics simulations of molecular solutions for system sizes of several thousands of atoms. As example of production applications we applied the method to molecular dynamics simulations in the isobaric-isothermal ensemble of the binary mixtures cyclohexane-methanol and toluene-methanol at different molar fractions of methanol. This record contains all CP2K input files necessary for the MD simulations, as well as 30ps trajectory files.