A rigorous theoretical treatment of vibrational energy relaxation in solution has been developed based on a general theory of dynamics of chemical reactions in solution. Algorithms which permit the construction of a physically realistic generalized Langevin equation of motion for the energy relaxation dynamics of a specified normal mode coordinate immersed at infinite dilution in monatomic and molecular solvent are developed. These algorithms permit the construction, from equilibrium solute-solvent pair correlation functions, of the liquid state frequency of the normal mode, \omega\sb{l}, and of the Gaussian model approximation to the autocorrelation function \langle\tilde{\cal F}(t)\tilde{\cal F}\rangle\sb{o} of the fluctuating force exerted by the solvent on the solute normal mode. From these quantities, one may compute the vibration energy relaxation time T\sb1 of the solute normal mode and assess the relative importance of the various energy dissipation pathways, solute vibration β solvent, solute vibration β solute translation β solvent, and solute vibration β solute translation, solute rotation β solvent. Numerical studies are presented for the prototype case of a diatomic solute in a monatomic solvent along with studies of more chemically interesting systems in molecular solvents. The study of VER in neat O\sb2 is presented