An efficient approach for an accurate quantum mechanical (QM) modeling of infrared (IR) spectra of condensed-phase systems is described. An ensemble of energetically low-lying cluster structures of a solute molecule surrounded by an explicit shell of solvent molecules is efficiently generated at the semiempirical tight-binding QM level and then reoptimized at the density functional theory level of theory. The IR spectrum of the solvated molecule is obtained as a thermodynamic average of harmonically computed QM spectra for all significantly populated cluster structures. The accuracy of such simulations in comparison to experimental data for some organic compounds and their solutions is shown to be the same or even better than the corresponding QM computations of the gas-phase IR spectrum for the isolated molecule.