This paper deals with equilibrium statistical mechanical aspects of the extension of the MTGLE approach for chemical reaction dynamics in liquids [S. A. Adelman, Adv. Chem. Phys. 53, 611 (1983)] so that it is applicable to molecular as well as monatomic solvents. This extension is necessary in order to conveniently treat energy exchange between the solute molecules and the solvent vibrational degrees of freedom. The analysis yields a separation of equilibrium solvent density fluctuations into translation-rotational (TR) and vibrational (V) contributions. The TR fluctuations are treated within the rigid solvent molecule model. The V fluctuations are treated within the harmonic oscillator approximation. The analysis is carried out in terms of a set of generalized solvent phase space coordinates S=pvvpww where v and w are, respectively, the V and TR coordinates and where pv and pw are the corresponding conjugate momenta. In this coordinate system, canonical ensemble distribution function for the solvent given that the solute is fixed at configuration point r0, denoted by fCA[S;r0], may be factorized as fCA[S;r0] =fCA[ pvv] fCA[ pw w;v0r0] where fCA[ pww;v0r0] and fCA[pvv] are respectively: the rigid solvent model approximation to fCA[S;r0]; a vibrational phase space probability distribution function which governs equilibrium mean field V fluctuations in the pure solvent. This vibrational probability distribution function may be written as fCA[ pvv] =Z−1 exp[−β(Tv+W(v))] where Z=pure solvent partition function, where Tv =pure solvent kinetic energy, and where W[v] is a vibrational potential of mean force. The above factorization of fCA[S;r0] holds if: (i) Coriolis coupling between the v and w coordinates is ignored; (ii) The solvent vibrational frequencies are sufficiently high. Within the harmonic approximation to W[v], fCA[ pvv] =[2πkβT]−rNs det ω2 exp [−β(Tv+ 1/2 yTω2y)] where y are the mass-weighted displacements of the v coordinates from equilibrium, where rNs =number of solvent normal modes, and where ω2 is an (rNs×rNs) -dimensional dynamical matrix which determines the vibrational frequency spectrum of the pure solvent within the harmonic approximation. This dynamical matrix may be constructed from the canonical ensemble distribution function of the pure rigid molecule solvent. Thus given the harmonic approximation to W[v] the canonical ensemble distribution function of the nonrigid solvent fCA[S;r0] may be constructed from the equilibrium properties of the rigid molecule solvent.