We adapt the three-dimensional reference interaction site model (3D-RISM) to calculate the potentials of mean force for ion-molecular solution as a difference between the chemical potential of solvation of a cluster of solutes and of individual ones. The method yields the solvation structure around the cluster of solutes in detail. The solvation chemical potential is obtained for the three-dimensional hypernetted chain (3D-HNC) closure as well as for its partial linearization (3D-PLHNC approximation). The solvation chemical potential is obtained in a closed analytical form for both the 3D-HNC and 3D-PLHNC closures. The 3D-RISM integral equations are solved by using the supercell technique. A straightforward supercell treatment of ionic solute in polar molecular solvent leads to a big error in the potential of mean force as well as the solvation chemical potential, which for simple ions in water amounts to about 35 kcal/mol. We elaborated corrections to the 3D-RISM integral equations, alleviating the artifact of the supercell periodicity with an accuracy of 0.05 kcal/mol or better and restoring the long-range asymptotics of the solute-solvent correlation functions. The dielectrically consistent site-site RISM/HNC theory (DRISM/HNC) is employed for the solvent correlations to provide a proper description of the dielectric properties of solution. This allowed us to extend the description to solution at a finite salt concentration. We converge both the 3D-RISM and site-site DRISM integral equations by using the method of modified direct inversion in the iterative subspace. Owing to the proper initial guess of the correlation functions, iteration begins at once for a given temperature and full molecular charge, avoiding a gradual decrease of the temperature and increase of the site charges, which greatly reduces the computation time. We calculate and discuss the potentials of mean force for sodium chloride in ambient water at infinite dilution as well as at a finite concentration.