The orientationally dependent potential of mean force (PMF) between two charged polyatomic solutes immersed in a polar molecular solvent is obtained by using the three-dimensional reference interaction site model (3D RISM) of the integral equation theory. We propose a partially linearized hypernetted chain (PLHNC) closure. The method is applied to the N,N-dimethylaniline cation (DMA(+)) and the anthracene anion (AN(-)) in acetonitrile solvent (CH3CN). We solve the 3D RISM integral equations by employing the modified direct inversion in the iterative subspace (MDIIS) method, The 3D site distributions of Solvent around each solute are obtained and discussed. The PMF between the solutes is calculated as a 3D profile dependent on the relative position of the solutes at six characteristic relative orientations. The PMF obtained is in qualitative agreement with the results of molecular dynamics simulations. In the solvent, the AN(-) solute effectively attracts the DMA(+) dimethylamino group and repels its phenyl ring. We found that the most stable relative arrangement of the DMA(+) and AN(-) molecules in acetonitrile solvent is different from that in gas phase.