1D-RISM

The one-dimensional reference interaction site model (1D-RISM) is an integral-equation theory that describes the structure of a homogeneous bulk solution through radial correlation functions between the atomic sites of its molecules. Because it averages over molecular orientations, every correlation function depends only on the site–site distance \(r\). In the ESM-RISM workflow, 1D-RISM supplies the solvent susceptibility \(\chi_{\alpha\gamma}\), which describes solvent–solvent correlations and is the input to 3D-RISM and the ESM-RISM method. We use Hartree atomic units (\(\hbar = e = m_\mathrm{e} = 1\)) and \(\beta = 1/k_\mathrm{B}T\) throughout.

Ornstein–Zernike equation for molecular liquids

For a simple liquid, the Ornstein–Zernike (OZ) equation relates the total correlation function \(h(r) = g(r) - 1\) to the direct correlation function \(c(r)\):

\[\begin{equation} h(r_{12}) = c(r_{12}) + \rho \int \mathrm{d}\boldsymbol{r}_3\, c(r_{13})\, h(r_{32}) \label{eq:oz} \end{equation}\]

Here \(\rho\) is the number density and \(g(r)\) the radial distribution function. The equation splits the correlation between particles 1 and 2 into a direct part \(c\) and a chain of indirect correlations mediated by a third particle. In Fourier space the convolution becomes a product, \(h(g) = c(g)/[1 - \rho c(g)]\).

For a molecular liquid, combining the position \(\boldsymbol{r}_i\) and orientation \(\Omega_i\) into \(i = (\boldsymbol{r}_i, \Omega_i)\) gives the molecular OZ equation:

\[\begin{equation} h(1,2) = c(1,2) + \frac{\rho}{8\pi^2} \int \mathrm{d}(3)\, c(1,3)\, h(3,2) \label{eq:moz} \end{equation}\]

where \(\mathrm{d}(3) = \mathrm{d}\boldsymbol{r}_3\,\mathrm{d}\Omega_3\) and \(\int \mathrm{d}\Omega = 8\pi^2\) for nonlinear molecules. The correlation functions now depend on six variables, which makes a direct numerical solution impractical.

Site–site representation: the RISM equation

Chandler and Andersen treated a molecule as a set of interaction sites (usually atoms) and described the liquid with orientationally averaged site–site correlation functions \(h_{\alpha\gamma}(r)\) and \(c_{\alpha\gamma}(r)\). The Greek indices \(\alpha, \gamma\) run over the sites of all solution components (solvent molecules and ions). Hirata and Rossky extended the theory to polar liquids made of charged sites (extended RISM, XRISM). In matrix form, the RISM equation reads

\[\begin{equation} \begin{split} h &= \omega * c * \omega + \omega * c * \rho h \\ &= [1 - \omega * c\rho]^{-1}\, \omega * c * \omega \end{split} \label{eq:rism} \end{equation}\]

where \(h\), \(c\), and \(\omega\) are matrices in the site indices, \(\rho\) is the diagonal matrix of number densities \(\rho_\alpha\) (the density of the component carrying site \(\alpha\)), and \(*\) denotes a convolution. We evaluate the convolutions as matrix products in Fourier space.

Intramolecular correlation function

The intramolecular correlation function \(\omega_{\alpha\gamma}\) encodes the molecular geometry. When sites \(\alpha\) and \(\gamma\) belong to the same molecule at a fixed distance \(l_{\alpha\gamma}\),

\[\begin{equation} \omega_{\alpha\gamma}(g) = \frac{\sin(g l_{\alpha\gamma})}{g l_{\alpha\gamma}} \label{eq:omega} \end{equation}\]

This is the Fourier transform of a spherical shell of radius \(l_{\alpha\gamma}\), \(\omega_{\alpha\gamma}(r) = \delta(r - l_{\alpha\gamma})/(4\pi l_{\alpha\gamma}^2)\); for \(\alpha = \gamma\) (\(l_{\alpha\alpha} = 0\)) it reduces to \(\omega_{\alpha\alpha}(g) = 1\). When \(\alpha\) and \(\gamma\) belong to different molecules, \(\omega_{\alpha\gamma} = 0\). A monatomic ion is a one-site molecule with only \(\omega_{\alpha\alpha} = 1\).

Solvent susceptibility

The right-hand side of the first line of Eq. \(\eqref{eq:rism}\) can be written as \(\omega * c * (\omega + \rho h)\). The quantity in parentheses is the solvent susceptibility:

\[\begin{equation} \chi_{\alpha\gamma}(g) = \omega_{\alpha\gamma}(g) + \rho_\alpha h_{\alpha\gamma}(g) \label{eq:chi} \end{equation}\]

\(\chi_{\alpha\gamma}\) describes how the density of site \(\alpha\) responds to a perturbation acting on site \(\gamma\), including both intramolecular (\(\omega\)) and intermolecular (\(\rho h\)) contributions. This function, built from the converged 1D-RISM solution, is what 3D-RISM and ESM-RISM use.

Interaction potential and closures

The RISM equation contains two unknown functions, \(h\) and \(c\), so a second relation, the closure, is needed. The exact closure is

\[\begin{equation} g_{\alpha\gamma}(r) = \exp\left[-\beta u_{\alpha\gamma}(r) + h_{\alpha\gamma}(r) - c_{\alpha\gamma}(r) + b_{\alpha\gamma}(r)\right] \label{eq:closure-exact} \end{equation}\]

where \(g_{\alpha\gamma} = h_{\alpha\gamma} + 1\) is the pair distribution function and \(b_{\alpha\gamma}\) is the bridge function, which contains higher-order correlations. Since the bridge function is generally unknown, we approximate it.

The site–site pair potential \(u_{\alpha\gamma}\) is the sum of a 12-6 Lennard-Jones (LJ) potential and the Coulomb potential:

\[\begin{equation} u_{\alpha\gamma}(r) = 4\varepsilon_{\alpha\gamma}\left[\left(\frac{\sigma_{\alpha\gamma}}{r}\right)^{12} - \left(\frac{\sigma_{\alpha\gamma}}{r}\right)^{6}\right] + \frac{q_\alpha q_\gamma}{r} \label{eq:pair} \end{equation}\]

where \(q_\alpha\) is the site charge. LJ parameters between unlike sites follow the Lorentz–Berthelot rules, \(\sigma_{\alpha\gamma} = (\sigma_\alpha + \sigma_\gamma)/2\) and \(\varepsilon_{\alpha\gamma} = \sqrt{\varepsilon_\alpha \varepsilon_\gamma}\).

HNC closure

Neglecting the bridge function gives the hypernetted-chain (HNC) closure:

\[\begin{equation} g_{\alpha\gamma}(r) = \exp\left[-\beta u_{\alpha\gamma}(r) + h_{\alpha\gamma}(r) - c_{\alpha\gamma}(r)\right] \label{eq:hnc} \end{equation}\]

HNC works well for systems with strong electrostatic interactions, but near strongly attractive sites the exponential grows too large and the iteration tends to diverge.

KH closure

Kovalenko and Hirata proposed a closure (the KH closure) that linearizes HNC wherever the exponent is positive:

\[\begin{equation} g_{\alpha\gamma}(r) = \begin{cases} \exp\left[d_{\alpha\gamma}(r)\right] & d_{\alpha\gamma}(r) \le 0 \\ 1 + d_{\alpha\gamma}(r) & d_{\alpha\gamma}(r) > 0 \end{cases}, \qquad d_{\alpha\gamma}(r) = -\beta u_{\alpha\gamma}(r) + h_{\alpha\gamma}(r) - c_{\alpha\gamma}(r) \label{eq:kh} \end{equation}\]

It coincides with HNC where the distribution is depleted (\(d \le 0\)) and grows only linearly where the distribution accumulates. This keeps the solution stable near metal surfaces and charged solutes. ESM-RISM calculations mainly use the KH closure.

Separating the long-range Coulomb interaction

At large distances the direct correlation function behaves as \(c_{\alpha\gamma}(r) \to -\beta u_{\alpha\gamma}(r)\), so it carries the \(1/r\) tail of the Coulomb interaction. Fourier-transforming this tail numerically degrades accuracy and stability, so we split \(c\) into short- and long-range parts:

\[\begin{equation} c_{\alpha\gamma}(r) = c^\mathrm{S}_{\alpha\gamma}(r) + c^\mathrm{L}_{\alpha\gamma}(r), \qquad c^\mathrm{L}_{\alpha\gamma}(r) = -\beta u^\mathrm{L}_{\alpha\gamma}(r) \label{eq:c-split} \end{equation}\]

The long-range part is the Coulomb potential screened by a Gaussian of width \(\tau\):

\[\begin{equation} u^\mathrm{L}_{\alpha\gamma}(r) = q_\alpha q_\gamma \frac{\operatorname{erf}(r/\tau)}{r}, \qquad u^\mathrm{L}_{\alpha\gamma}(g) = \frac{4\pi q_\alpha q_\gamma}{g^2} \exp\left(-\frac{g^2\tau^2}{4}\right) \label{eq:uL} \end{equation}\]

Because \(u^\mathrm{L}\) is analytic, we compute it once before the iteration and optimize only the short-range part \(c^\mathrm{S}\). With \(u^\mathrm{S} = u - u^\mathrm{L}\), the identity \(-\beta u + h - c = -\beta u^\mathrm{S} + h - c^\mathrm{S}\) holds, so the closures \(\eqref{eq:hnc}\) and \(\eqref{eq:kh}\) involve only short-range quantities. We solve the RISM equation \(\eqref{eq:rism}\) in Fourier space after adding the analytic Fourier components of \(c^\mathrm{L}\). A typical choice is \(\tau \approx 1\) Å.

Dielectrically consistent RISM (DRISM)

XRISM systematically underestimates the relative permittivity \(\varepsilon\) of solutions of dipolar molecules. The deficiency stems from neglecting the bridge function, which carries higher-order correlations such as dipole–dipole correlations. The dielectrically consistent RISM (DRISM) of Perkyns and Pettitt takes \(\varepsilon\) as an empirical parameter and modifies the RISM equation so that Coulomb interactions between charges in the solution are screened correctly:

\[\begin{equation} h' = [1 - \omega' * c\rho]^{-1}\, \omega' * c * \omega', \qquad h' = h - \xi, \quad \omega' = \omega + \rho\xi \label{eq:drism} \end{equation}\]

The dielectric correction function \(\xi_{\alpha\gamma}\) is nonzero only when both \(\alpha\) and \(\gamma\) are sites of the polar molecule; it vanishes for any pair that involves an ion:

\[\begin{equation} \xi_{\alpha\gamma}(g) = j_0(g d^x_\alpha)\, j_0(g d^y_\alpha)\, j_1(g d^z_\alpha)\; \tilde{h}(g)\; j_0(g d^x_\gamma)\, j_0(g d^y_\gamma)\, j_1(g d^z_\gamma) \label{eq:xi} \end{equation}\]

Here \(j_0\) and \(j_1\) are spherical Bessel functions of order zero and one, and \(d^i_\alpha\) (\(i = x, y, z\)) are the coordinates of site \(\alpha\) measured from the center of absolute charges of the polar molecule, oriented with its dipole along \(z\). The kernel \(\tilde{h}(g)\) is

\[\begin{equation} \tilde{h}(g) = \frac{1}{\rho_\mathrm{tot}}\left(\frac{\varepsilon - 1}{y} - 3\right) \exp\left(-\frac{a^2 g^2}{4}\right), \qquad y = \frac{4\pi \rho_\mathrm{pol} \mu^2 \beta}{9} \label{eq:htilde} \end{equation}\]

where \(\rho_\mathrm{tot}\) is the total number density, \(\rho_\mathrm{pol}\) the number density of the polar molecule, \(\mu\) its dipole moment, and \(y\) the dipole density. The parameter \(a\) restricts the correction to long wavelengths; about 0.5 Å is used for aqueous solutions. In the uncorrelated limit \(\varepsilon - 1 = 3y\), the correction vanishes.

Because \(\omega' + \rho h' = \omega + \rho h\), Eq. \(\eqref{eq:chi}\) still yields \(\chi_{\alpha\gamma}\) from the DRISM solution. The dielectric correction therefore passes directly into 3D-RISM and ESM-RISM through the solvent susceptibility. Hagiwara et al. combined DRISM with ESM-RISM and showed that it improves on XRISM for the mean activity coefficients of aqueous NaCl, and for the potential of zero charge and the double-layer capacitance of the Pt(111)/aqueous HCl interface.

Excess chemical potential

From the 1D-RISM solution, the excess chemical potential (solvation free energy) of a component \(s\) follows in closed form. For the KH closure,

\[\begin{equation} \Delta\mu_s = k_\mathrm{B}T \sum_{\alpha \in s} \sum_\gamma 4\pi\rho_\gamma \int_0^\infty r^2\, \mathrm{d}r \left[\frac{1}{2} h_{\alpha\gamma}^2(r)\, \Theta\bigl(-h_{\alpha\gamma}(r)\bigr) - c_{\alpha\gamma}(r) - \frac{1}{2} h_{\alpha\gamma}(r)\, c_{\alpha\gamma}(r)\right] \label{eq:mu-kh} \end{equation}\]

where \(\Theta\) is the Heaviside step function. Removing \(\Theta\) gives the HNC (Singer–Chandler) expression, and dropping the \(\frac{1}{2}h^2\) term gives the Gaussian fluctuation (GF) expression. These expressions are consistent with their closures and yield the free energy from a single calculation without thermodynamic integration, a key advantage of RISM theory.

For electrolyte solutions at finite concentration, the mean activity coefficient \(\gamma_{m,\pm}\) follows from the excess chemical potentials:

\[\begin{equation} \nu k_\mathrm{B}T \ln \gamma_{m,\pm} = k_\mathrm{B}T \ln \frac{\rho_\mathrm{w}}{\rho_{\mathrm{w},\infty}} + \sum_i s_i\, \Delta\mu_{\mathrm{ex},i} \label{eq:activity} \end{equation}\]

Here \(i\) labels the ionic species, \(s_i\) are stoichiometric coefficients with \(\nu = \sum_i s_i\), and \(\rho_\mathrm{w}\) and \(\rho_{\mathrm{w},\infty}\) are the number densities of water at finite concentration and at infinite dilution. XRISM strongly underestimates \(\gamma_{m,\pm}\), whereas DRISM reproduces experiment well.

Numerical solution

Since the 1D-RISM correlation functions are radial, the Fourier transform reduces to a spherical Bessel (sine) transform:

\[\begin{equation} f(g) = \frac{4\pi}{g}\int_0^\infty r \sin(gr)\, f(r)\, \mathrm{d}r, \qquad f(r) = \frac{1}{2\pi^2 r}\int_0^\infty g \sin(gr)\, f(g)\, \mathrm{d}g \label{eq:hankel} \end{equation}\]

The calculation iterates as follows:

1. Start from an initial \(c^\mathrm{S}_{\alpha\gamma}(r)\), Fourier-transform it, and add \(c^\mathrm{L}\).
2. Solve the RISM equation \(\eqref{eq:rism}\) (or the DRISM equation \(\eqref{eq:drism}\)) in Fourier space as a matrix operation at each \(g\) to obtain \(h_{\alpha\gamma}(g)\).
3. Back-transform to \(h_{\alpha\gamma}(r)\) and insert it into the closure to obtain a new \(c^\mathrm{S}_{\alpha\gamma}(r)\).
4. Repeat until the residual falls below the convergence threshold.

Simple mixing converges slowly, so we accelerate the iteration with the modified direct inversion in the iterative subspace (MDIIS) method. Hagiwara et al. used a real-space cutoff of 999.8 Bohr, a reciprocal-space cutoff of 18.763 Bohr\(^{-1}\), 5973 grid points, and an MDIIS convergence threshold of \(10^{-8}\).

Passing the result to 3D-RISM and ESM-RISM

The main output of 1D-RISM is the solvent susceptibility \(\chi_{\alpha\gamma}(g)\), assembled from the converged \(h_{\alpha\gamma}(g)\) through Eq. \(\eqref{eq:chi}\). It depends only on the bulk solution, not on the solute.

• 3D-RISM interpolates \(\chi_{\alpha\gamma}(|\boldsymbol{g}|)\) onto three-dimensional reciprocal lattice vectors and solves for the three-dimensional distribution around the solute.
• The ESM-RISM method inverse-Fourier-transforms \(\chi_{\alpha\gamma}(g)\) along the surface normal \(g_z\) to obtain the Laue representation \(\chi_{\alpha\gamma}(g_\parallel, |z - z'|)\). Because \(\omega_{\alpha\gamma}(g)\) does not decay at large \(g\), this transform needs special care (see the ESM-RISM chapter).

The solvent model (site geometry, LJ parameters, charges), temperature, concentration, and the choice of dielectric correction are all fixed at this stage, so the accuracy of 1D-RISM governs the results of 3D-RISM and ESM-RISM.

References

• D. Chandler and H. C. Andersen, "Optimized cluster expansions for classical fluids. II. Theory of molecular liquids," J. Chem. Phys. 57, 1930 (1972). doi
• F. Hirata and P. J. Rossky, "An extended RISM equation for molecular polar fluids," Chem. Phys. Lett. 83, 329 (1981). doi
• F. Hirata, B. M. Pettitt, and P. J. Rossky, "Application of an extended RISM equation to dipolar and quadrupolar fluids," J. Chem. Phys. 77, 509 (1982). doi
• T. Morita, "Theory of classical fluids: Hyper-netted chain approximation, I," Prog. Theor. Phys. 20, 920 (1958). doi
• A. Kovalenko and F. Hirata, "Self-consistent description of a metal–water interface by the Kohn–Sham density functional theory and the three-dimensional reference interaction site model," J. Chem. Phys. 110, 10095 (1999). doi
• J. S. Perkyns and B. M. Pettitt, "A dielectrically consistent interaction site theory for solvent–electrolyte mixtures," Chem. Phys. Lett. 190, 626 (1992). doi
• J. Perkyns and B. M. Pettitt, "A site–site theory for finite concentration saline solutions," J. Chem. Phys. 97, 7656 (1992). doi
• I. S. Joung, T. Luchko, and D. A. Case, "Simple electrolyte solutions: Comparison of DRISM and molecular dynamics results for alkali halide solutions," J. Chem. Phys. 138, 044103 (2013). doi
• S. J. Singer and D. Chandler, "Free energy functions in the extended RISM approximation," Mol. Phys. 55, 621 (1985). doi
• R. M. Levy, M. Belhadj, and D. B. Kitchen, "Gaussian fluctuation formula for electrostatic free-energy changes in solution," J. Chem. Phys. 95, 3627 (1991). doi
• A. Kovalenko, S. Ten-no, and F. Hirata, "Solution of three-dimensional reference interaction site model and hypernetted chain equations for simple point charge water by modified method of direct inversion in iterative subspace," J. Comput. Chem. 20, 928 (1999). doi
• S. Nishihara and M. Otani, "Hybrid solvation models for bulk, interface, and membrane: Reference interaction site methods coupled with density functional theory," Phys. Rev. B 96, 115429 (2017). doi
• S. Hagiwara, S. Nishihara, F. Kuroda, and M. Otani, "Development of a dielectrically consistent reference interaction site model combined with the density functional theory for electrochemical interface simulations," Phys. Rev. Mater. 6, 093802 (2022). doi