We have developed a self-consistent description of a metal - molecular liquid interface by combination of the Kohn-Sham density functional theory (KS DFT) for the electronic density, and the three-dimensional reference interactions site model (3D RISM) integral equation theory for the classical site distribution profiles of molecular liquid. The electron and classical subsystems are coupled in the mean field approximation. The classical potentials of the metal acting on species of the liquid are taken in the linear response regime. Many-body effects of dense liquid on metal valence electrons are allowed for by averaging pseudopotentials of liquid molecules over the site distributions of liquid. The coupled KS DFT and 3D RISM equations are solved simultaneously by using the procedure of dynamic relaxation. The proposed approach is substantially less time-consuming as compared to a Car-Parrinello type simulation. A partial linearization of the hypernetted chain (PLHNC) closure to the RISM equation is proposed. The calculation is performed in the supercell technique for water at normal conditions in contact with the (100) FCC slab of a metal roughly modeled after copper. The results are in good agreement with the Car-Parrinello simulation for the same model. Further applications of the method proposed are discussed.