The derivation of a three-dimensional integral equation for solute molecule-solvent site correlation functions is presented. The equation is obtained by averaging the Ornstein-Zernicke equation for molecular liquids over orientations of the solvent molecule consistent with one site of the solvent remaining at a fixed distance from a solute-based origin. The approach is similar to that adopted in the reduction leading to the reference interaction site model (RISM) equations but retains full three-dimensional information regarding the structure of the reference solute molecule. The proposed equation can be solved using three-dimensional HNC-like closures, of which three different forms are discussed. A formulation which allows the introduction of long range interactions through a renormalization of the equation is also presented. Applications to various molecular liquids indicate that the proposed theory provides pair correlation functions that are in better agreement with molecular dynamics simulations than those obtained using the extended RISM formulation. Furthermore, qualitative errors in the correlation functions, frequently seen in results from RISM calculations are completely eliminated through geometrical averaging of the Mayer function in the 3D HNC closure. Prospects for the development of a novel mean field theory of solvation are also discussed. (C) 1997 American Institute of Physics.