Accurate descriptions of intermolecular interactions are of great importance in simulations of molecular liquids. We present an electronic structure method that combines the accuracy of the Harris functional approach with the computational efficiency of approximately linear-scaling density functional theory (DFT). The proposed method allows for simulations with accuracies close to the Kohn-Sham DFT reference. 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, and benchmark data of the energy correction.