The three-dimensional reference interaction site model with the closure relation by Kovalenko and Hirata (3D-RISM-KH) in combination with the density functional theory (DFT) method has been implemented in the Amsterdam density functional (ADF) software package. The analytical first derivatives of the free energy with respect to displacements of the solute nuclear coordinates have also been developed. This enables study of chemical reactions, including reaction coordinates and transition state search, with the molecular solvation described from the first principles. The method yields all of the features available by using other solvation approaches, for instance infrared spectra of solvated molecules. To evaluate the accuracy of the present method, test calculations have been carried out for a number of small molecules, including four glycine conformers, a set of small organic compounds, and carbon nanotubes of various lengths in aqueous solution. Our predictions for the solvation free energy agree well with other approaches as well as experiment. This new development makes it possible to calculate at modest computational cost the electronic properties and molecular solvation structure of a solute molecule in a given molecular liquid or mixture from the first principles.