In this paper, we present a combined approach of quantum mechanics and statistical mechanical integral equation theory in three dimensions to describe the solute-solvent electrostatic interaction and mutual polarization in an infinitely dilute solution of a polar solute molecule in polar solvent. In this combined method, the integral equation theory 3d-RISM-HNC (J. Phys. Chem. B 2000, 104, 796) provided the average solvent particle densities and solvent charge distribution surrounding a solute molecule. The Schrodinger equation of the solute molecule in solvent background charges was solved with the quantum mechanical SCF Hartree-Fock (HF) method. The electrostatic potentials were calculated using the electronic density of the solute, which can be fitted to the atomic ESP charges. Those charges were, in turn, used in 3d-RISM-HNC calculations, which were solved iteratively. We have seen a significant increase in the solute dipole moment and the solute-solvent electrostatic interaction energy. In the case study of N-methyl amine in liquid water, we observed higher hydrogen density of water near O(CO) and higher oxygen density of water near H(NH). It has been shown that the mutual polarization between solute and solvent has significant effect on the thermodynamic properties.