A set of integral equations is derived for site-site pair correlation functions for diatomic fluids from the Kirkwood hierarchy of correlation functions. The integral equations obtained are self-consistent field equations which are generalizations of the Percus-Yevick integral equation for simple liquids. In this theory the intra- and intermolecular correlation functions are self-consistently coupled. The integral equations reduce to the integral equations in the reference site interaction model (RISM) theory of diatomic liquids if an approximation is made to the direct correlation functions and the molecule is assumed to be rigid. The theory is also generalized to polyatomic fluids. The predictions of the theory are compared with those by the RISM theory and Monte Carlo simulations available in the literature. The computational requirement for the integral equations of the present theory is comparable to that of the RISM theory, but its predictions are more accurate than those by RISM theory at low densities and are practically the same as the latter at high densities.