We developed and implemented a statistical mechanical integral equation theory to describe the hydration structure of complex molecules. The theory, which is an extension of the reference interaction site model (RISM) in three dimensions, yields the average density from the solvent interactions sites at all points r around a molecular solute of arbitrary shape. Both solute-solvent electrostatic and van der Waals interactions are fully included, and solvent packing is taken into account. The approach is illustrated by calculating the average oxygen and hydrogen density of liquid water around two molecular solutes : water and N-methylacetamide. Molecular dynamics simulations are performed to test the results obtained from the integral equation. It is observed that important microscopic structural features of the average water density due to hydrogen bonding are reproduced by the integral equation. The integral equation has a simple formal structure and is easy to implement numerically. It offers a powerful alternative to computer simulations with explicit solvent molecules and to continuum solvent representations for incorporating solvation effects in a wide range of applications.