A new approach to investigate a molecular recognition process of protein is presented based on the three-dimensional reference interaction site model (3D-RISM) theory, a statistical mechanics theory of molecular liquids. Numerical procedure for solving the conventional 3D-RISM equation consists of two steps. In step 1, we solve ordinary RISM (or 1D-RISM) equations for a solvent mixture including target ligands in order to obtain the density pair correlation functions (PCF) among molecules in the solution. Then, we solve the 3D-RISM equation for a solute-solvent system to find three-dimensional density distribution functions (3D-DDF) of solvent species around a protein, using PCF obtained in the first step. A key to the success of the method was to regard a target ligand as one of "solvent" species. However, the success is limited due to a difficulty of solving the 1D-RISM equation for a solvent mixture, including large ligand molecules. In the present paper, we propose a method which eases the limitation concerning solute size in the conventional method. In this approach, we solve a solute-solute 3D-RISM equations for a protein-ligand system in which both proteins and ligands are regarded as "solutes" at infinite dilution. The 3D- and 1D-RISM equations are solved for protein-solvent and ligand-solvent systems, respectively, in order to obtain the 3D- and 1D-DDF of solvent around the solutes, which are required for solving the solute-solute 3D-RISM equation. The method is applied to two practical and noteworthy examples concerning pharmaceutical design. One is an odorant binding protein in the Drosophila melanogaster , which binds an ethanol molecule. The other is phospholipase A2, which is known as a receptor of acetylsalicylic acid or aspirin. The result indicates that the method successfully reproduces the binding mode of the ligand molecules in the binding sites measured by the experiments.