The molecules of real liquids have internal degrees of freedom that may couple with the external coordinates of position and orientation so that they affect and are affected by the microscopic liquid structure, For cases where the internal coordinates possess a Boltzmann-like distribution, a procedure was recently proposed [Phys. Rev. E 55, 426 (1997)] whereby the internal coordinates are incorporated into the conventional integral equation formulation of classical liquid state theory with no approximations beyond some reliable closure relation familiar from simple liquids. The basis of the procedure is expansions in special orthogonal polynomials of the internal coordinates. Here we use this technique to obtain the structural, thermodynamic, and electrostatic properties of a classical liquid of polar polarizable molecules, with classical Drude oscillators modeling the infernal variable of fluctuating polarization, Sample results obtained using several approximate closures are compared with simulation data.