A combined method of the Dirac–Hartree–Fock (DHF) method and the reference interaction-site model (RISM) theory is reported; this is the initial implementation of the coupling of the four-component relativistic electronic structure theory and an integral equation theory of molecular liquids. In the method, the DHF and RISM equations are solved self-consistently, and therefore the electronic structure of the solute, including relativistic effects, and the solvation structure are determined simultaneously. The formulation is constructed based on the variational principle with respect to the Helmholtz energy, and analytic free energy gradients are also derived using the variational property. The method is applied to the iodine ion (I), methyl iodide (CH3I), and hydrogen chalcogenide (H2X, where X = O–Po) in aqueous solutions, and the electronic structures of the solutes, as well as the solvation free energies and their component analysis, solvent distributions, and solute–solvent interactions, are discussed.