3D-RISM

The three-dimensional reference interaction site model (3D-RISM) is an integral-equation theory that yields the distribution of solvent sites around a solute as three-dimensional functions. Whereas 1D-RISM treats the solute through orientationally averaged radial functions, 3D-RISM fixes the solute in space and retains its full shape and charge distribution. When the solute is an electron–ion system described by density functional theory (DFT), 3D-RISM becomes a hybrid solvation model that couples DFT and liquid-state theory self-consistently. We use Hartree atomic units (\(\hbar = e = m_\mathrm{e} = 1\)) and \(\beta = 1/k_\mathrm{B}T\), and denote solvent sites by \(\alpha, \gamma\) and solute atoms by \(I\).

From 1D-RISM to 3D-RISM

Consider a mixture of a solute (indices \(a, b\)) and a solvent (indices \(\alpha, \gamma, \lambda\)), with the solute at infinite dilution. The solute–solvent block of the site–site RISM equation \(h = \omega * c * \chi\) reads

\[\begin{equation} h_{a\gamma}(\boldsymbol{r}) = \sum_{\mu, \lambda} \int \mathrm{d}\boldsymbol{r}' \int \mathrm{d}\boldsymbol{r}''\, \omega_{a\mu}(\boldsymbol{r} - \boldsymbol{r}')\, c_{\mu\lambda}(\boldsymbol{r}' - \boldsymbol{r}'')\, \chi_{\lambda\gamma}(\boldsymbol{r}'') \label{eq:rism-uv} \end{equation}\]

where \(\mu, \lambda\) run over all solute and solvent sites. Instead of representing the solute structure through an orientational average with intramolecular correlation functions, 3D-RISM fixes the solute in space and stores its structure directly in three-dimensional correlation functions. We therefore set, for solute sites,

\[\begin{equation} \omega_{a\mu}(\boldsymbol{r}) = \delta_{a\mu}\, \delta(\boldsymbol{r}) \qquad (a \in \text{solute}) \label{eq:omega-solute} \end{equation}\]

and set the direct correlation between solute sites to zero. The sum over \(\mu\) then reduces to \(\mu = a\), and the sum over \(\lambda\) runs only over solvent sites. Treating the whole solute as a single three-dimensional object, we drop the index \(a\) and approximate the solvent–solvent susceptibility \(\chi_{\alpha\gamma}\) by the isotropic bulk function from 1D-RISM, unperturbed by the solute. This yields the 3D-RISM equation:

\[\begin{equation} h_\gamma(\boldsymbol{r}) = \sum_\alpha \int \mathrm{d}\boldsymbol{r}'\, c_\alpha(\boldsymbol{r}')\, \chi_{\alpha\gamma}(|\boldsymbol{r} - \boldsymbol{r}'|) \label{eq:3drism} \end{equation}\]

Here \(h_\gamma(\boldsymbol{r}) = g_\gamma(\boldsymbol{r}) - 1\) and \(c_\gamma(\boldsymbol{r})\) are the three-dimensional total and direct correlation functions of solvent site \(\gamma\) around the solute, and \(g_\gamma(\boldsymbol{r})\) is the density of site \(\gamma\) at \(\boldsymbol{r}\) normalized by its bulk density \(\rho_\gamma\). \(\chi_{\alpha\gamma}(g) = \omega_{\alpha\gamma}(g) + \rho_\alpha h_{\alpha\gamma}(g)\) is the solvent susceptibility from 1D-RISM. In Fourier space the convolution becomes a product:

\[\begin{equation} h_\gamma(\boldsymbol{g}) = \sum_\alpha c_\alpha(\boldsymbol{g})\, \chi_{\alpha\gamma}(|\boldsymbol{g}|) \label{eq:3drism-g} \end{equation}\]

Solute information enters only through \(c_\alpha\) and the interaction potential \(u_\gamma\) introduced below; solvent information enters only through \(\chi_{\alpha\gamma}\). Thanks to this separation, a single solvent calculation (1D-RISM) serves any number of solutes.

Solute–solvent interaction potential

The potential \(u_\gamma(\boldsymbol{r})\) acting on solvent site \(\gamma\) at \(\boldsymbol{r}\) is the sum of Lennard-Jones (LJ) interactions with the solute atoms and the electrostatic potential of the solute charge distribution:

\[\begin{equation} u_\gamma(\boldsymbol{r}) = \sum_I 4\varepsilon_{I\gamma}\left[\left(\frac{\sigma_{I\gamma}}{|\boldsymbol{r} - \boldsymbol{R}_I|}\right)^{12} - \left(\frac{\sigma_{I\gamma}}{|\boldsymbol{r} - \boldsymbol{R}_I|}\right)^{6}\right] - q_\gamma\, v_\mathrm{DFT}(\boldsymbol{r}) \label{eq:u3d} \end{equation}\]

Here \(\boldsymbol{R}_I\) is the position of solute atom \(I\) and \(q_\gamma\) the charge of the solvent site. LJ parameters are combined with the Lorentz–Berthelot rules. \(v_\mathrm{DFT}\) is the DFT electrostatic potential, the sum of the Hartree potential and the local part of the pseudopotentials:

\[\begin{equation} v_\mathrm{DFT}(\boldsymbol{r}) = v_\mathrm{H}(\boldsymbol{r}) + v_\mathrm{ion}(\boldsymbol{r}), \qquad v_\mathrm{H}(\boldsymbol{r}) = \int \mathrm{d}\boldsymbol{r}'\, \frac{n(\boldsymbol{r}')}{|\boldsymbol{r} - \boldsymbol{r}'|} \label{eq:vdft} \end{equation}\]

where \(n(\boldsymbol{r})\) is the electron density. Since \(v_\mathrm{DFT}\) is the potential energy of an electron (charge \(-1\)), a site of charge \(q_\gamma\) feels the electrostatic energy \(-q_\gamma v_\mathrm{DFT}\). A change in the electron distribution changes \(u_\gamma\) and hence the solvent distribution; this is how DFT couples to the solvent.

Closures

The 3D-RISM equation is solved together with a closure. The HNC closure reads

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

For solutes that strongly attract the solvent, such as metal electrodes, the HNC exponential becomes too large. We therefore use the KH closure, which Kovalenko and Hirata introduced for the metal–water interface:

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

Inside the solute, LJ repulsion drives \(g_\gamma \to 0\), so the solvent occupies only the region outside the solute.

Long-range asymptotics and charge handling

Far from the solute, \(c_\gamma(\boldsymbol{r}) \to -\beta u_\gamma(\boldsymbol{r})\), so a charged or dipolar solute gives \(c_\gamma\) a Coulomb-like tail. As in 1D-RISM, we split \(c_\gamma\) into short- and long-range parts:

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

The long-range part is the DFT electrostatic potential smeared by a Gaussian of width \(\tau\), given in Fourier space:

\[\begin{equation} u^\mathrm{L}_\gamma(\boldsymbol{g}) = -q_\gamma\, v_\mathrm{DFT}(\boldsymbol{g}) \exp\left(-\frac{g^2\tau^2}{4}\right) \label{eq:uL3d} \end{equation}\]

A typical choice is \(\tau \approx 1\) Å. For a given electron density, \(u^\mathrm{L}_\gamma\) is fixed before the 3D-RISM iteration, which then optimizes only the short-range part \(c^\mathrm{S}_\gamma\). The closure exponent involves only short-range quantities, \(-\beta u^\mathrm{S}_\gamma + h_\gamma - c^\mathrm{S}_\gamma\) with \(u^\mathrm{S}_\gamma = u_\gamma - u^\mathrm{L}_\gamma\), and the 3D-RISM equation \(\eqref{eq:3drism-g}\) uses \(c^\mathrm{S}_\alpha + c^\mathrm{L}_\alpha\).

For a charged solute, the total correlation function \(h_\gamma\) also has a long, Debye-screened tail. For isolated solutes, methods have been developed that separate the asymptotics of both \(h\) and \(c\) analytically. In a periodic supercell the \(\boldsymbol{g} = 0\) Coulomb component is undefined; as in ordinary plane-wave DFT, we set the average electrostatic potential to zero, and a net charge must be neutralized by a uniform background charge.

Solvation free energy

The 3D-RISM solution gives the solvation free energy (excess chemical potential) \(\Delta\mu_\mathrm{solv}\) in closed form. For the KH closure,

\[\begin{equation} \Delta\mu_\mathrm{solv} = k_\mathrm{B}T \sum_\gamma \rho_\gamma \int \mathrm{d}\boldsymbol{r} \left[\frac{1}{2} h_\gamma^2(\boldsymbol{r})\, \Theta\bigl(-h_\gamma(\boldsymbol{r})\bigr) - c_\gamma(\boldsymbol{r}) - \frac{1}{2} h_\gamma(\boldsymbol{r})\, c_\gamma(\boldsymbol{r})\right] \label{eq:mu-kh3d} \end{equation}\]

For the HNC closure, \(\Theta\) is removed. The Gaussian fluctuation (GF) approximation drops the \(\frac{1}{2}h^2\) term:

\[\begin{equation} \Delta\mu^\mathrm{GF}_\mathrm{solv} = k_\mathrm{B}T \sum_\gamma \rho_\gamma \int \mathrm{d}\boldsymbol{r} \left[- c_\gamma(\boldsymbol{r}) - \frac{1}{2} h_\gamma(\boldsymbol{r})\, c_\gamma(\boldsymbol{r})\right] \label{eq:mu-gf3d} \end{equation}\]

Hagiwara et al. obtained the distributions with the KH closure and evaluated the free energy with the GF expression. The KH and HNC expressions are stationary with respect to \(h_\gamma\) and \(c_\gamma\) at the solution; the next section uses this property to obtain the functional derivative with respect to the electron density.

Coupling to DFT

Helmholtz free energy

In DFT-RISM, the DFT total energy \(E_\mathrm{DFT}\) is replaced by the Helmholtz free energy \(A\), which adds the solvation free energy:

\[\begin{equation} A[n, \{\rho_\gamma\}] = E_\mathrm{DFT}[n] + \Delta\mu_\mathrm{solv}[n, \{\rho_\gamma\}] \label{eq:helmholtz} \end{equation}\]

\(E_\mathrm{DFT}[n]\) is the sum of the noninteracting kinetic energy, the external-potential energy, the Hartree energy, and the exchange-correlation energy. Minimizing \(A\) with respect to the electron density at fixed atomic positions and electron number gives the electronic structure of the solute in solution.

Feedback into the Kohn–Sham potential

The Kohn–Sham effective potential acquires a solvent contribution \(v_\mathrm{solv}\):

\[\begin{equation} v_\mathrm{eff}(\boldsymbol{r}) = v_\mathrm{H}(\boldsymbol{r}) + v_\mathrm{ion}(\boldsymbol{r}) + v_\mathrm{xc}(\boldsymbol{r}) + v_\mathrm{solv}(\boldsymbol{r}) \label{eq:veff} \end{equation}\]

Because the free energy \(\eqref{eq:mu-kh3d}\) is stationary with respect to \(h_\gamma\) and \(c_\gamma\), \(\delta\Delta\mu_\mathrm{solv}/\delta u_\gamma(\boldsymbol{r}) = \rho_\gamma g_\gamma(\boldsymbol{r})\). Equations \(\eqref{eq:u3d}\) and \(\eqref{eq:vdft}\) give \(\delta u_\gamma(\boldsymbol{r}')/\delta n(\boldsymbol{r}) = -q_\gamma/|\boldsymbol{r} - \boldsymbol{r}'|\), so

\[\begin{equation} v_\mathrm{solv}(\boldsymbol{r}) = \frac{\delta \Delta\mu_\mathrm{solv}}{\delta n(\boldsymbol{r})} = -\int \mathrm{d}\boldsymbol{r}'\, \frac{\rho_\mathrm{solv}(\boldsymbol{r}')}{|\boldsymbol{r} - \boldsymbol{r}'|}, \qquad \rho_\mathrm{solv}(\boldsymbol{r}) = \sum_\gamma q_\gamma \rho_\gamma g_\gamma(\boldsymbol{r}) \label{eq:vsolv} \end{equation}\]

\(\rho_\mathrm{solv}\) is the solvent charge density, and \(v_\mathrm{solv}\) is its electrostatic potential expressed as the potential energy of an electron. The LJ interaction does not depend on the electron density and therefore does not appear in \(v_\mathrm{solv}\). The calculation proceeds in a self-consistent loop:

1. Compute \(v_\mathrm{DFT}\) from the electron density \(n\) and build \(u_\gamma\) with Eq. \(\eqref{eq:u3d}\).
2. Solve the 3D-RISM equation with the closure to obtain \(g_\gamma\).
3. Compute \(\rho_\mathrm{solv}\) and \(v_\mathrm{solv}\) from Eq. \(\eqref{eq:vsolv}\), and solve the Kohn–Sham equations for a new \(n\).
4. Repeat until both \(n\) and \(g_\gamma\) converge.

Janak's theorem

With occupation numbers \(f_i\) and eigenvalues \(\epsilon_i\) of the Kohn–Sham orbitals, Janak's theorem also holds for the Helmholtz free energy \(\eqref{eq:helmholtz}\):

\[\begin{equation} \frac{\partial A[n, \{\rho_\gamma\}]}{\partial f_i} = \epsilon_i \label{eq:janak} \end{equation}\]

The proof follows the standard DFT one: \(\partial \Delta\mu_\mathrm{solv}/\partial f_i = \int \mathrm{d}\boldsymbol{r}\, |\psi_i(\boldsymbol{r})|^2 v_\mathrm{solv}(\boldsymbol{r})\) is absorbed as part of \(v_\mathrm{eff}\). The electron chemical potential is therefore \(\mu_\mathrm{e} = \partial A/\partial N\), which justifies evaluating the electrode potential in solution from the electron chemical potential.

Solving in a periodic supercell

When combined with plane-wave DFT, the 3D-RISM equation is solved in the same periodic supercell as DFT.

• Correlation functions are transformed by FFT between the real-space grid and the reciprocal lattice vectors \(\boldsymbol{g}\). Hagiwara et al. used a cutoff energy of 160 Ry for the correlation functions.
• \(\chi_{\alpha\gamma}(g)\), computed by 1D-RISM on a radial grid, is interpolated to the length \(|\boldsymbol{g}|\) of each reciprocal lattice vector and used in Eq. \(\eqref{eq:3drism-g}\).
• The closure is evaluated at each real-space grid point and the 3D-RISM equation at each \(\boldsymbol{g}\); the MDIIS method accelerates the iteration for \(c^\mathrm{S}_\gamma\).

This formulation is canonical. The number of solvent sites \(\gamma\) in the cell,

\[\begin{equation} N_\gamma = \rho_\gamma \int_\mathrm{cell} \mathrm{d}\boldsymbol{r}\, g_\gamma(\boldsymbol{r}) \label{eq:number} \end{equation}\]

remains constant during the calculation. This property is useful, for example, for evaluating the number of water molecules intercalated in layered materials.

Application to interfaces and its limitations

3D-RISM works well for molecules and ions in solution, but a periodic supercell makes it difficult to treat electrode/electrolyte interfaces.

• The slab and the solution between slabs repeat along the surface normal, so the solution region has a finite thickness and is not connected to a bulk reservoir. The inner potential of the bulk solution, which serves as the potential reference, cannot be defined.
• Under the canonical condition the number of solvent particles is fixed, so the theory cannot describe the formation of an electric double layer in which counter-ions enter or leave the solution as the electrode charges.
• A charged electrode requires a uniform background charge, which introduces artifacts into the interfacial electrostatic potential and energy.

The ESM-RISM method removes these problems. It rewrites the 3D-RISM equation in the Laue representation, Fourier-transformed only in the in-plane directions,

\[\begin{equation} h_\gamma(\boldsymbol{g}_\parallel, z) = \sum_\alpha \int \mathrm{d}z'\, c_\alpha(\boldsymbol{g}_\parallel, z')\, \chi_{\alpha\gamma}(\boldsymbol{g}_\parallel, |z - z'|) \label{eq:laue} \end{equation}\]

and imposes a boundary condition under which the solution extends to the bulk at infinity along the surface normal. The electrostatic potential is solved under open boundary conditions with the ESM method. Electrons and solvent particles can then both be treated in the grand-canonical ensemble, and the grand potential

\[\begin{equation} \Omega = A - \mu_\mathrm{e}(N - N_0) = A - \mu_\mathrm{e}\Delta N \label{eq:grand} \end{equation}\]

enables calculations at controlled electrode potential without a background charge. Here \(N_0\) is the electron number of the neutral system.

References

• D. Beglov and B. Roux, "An integral equation to describe the solvation of polar molecules in liquid water," J. Phys. Chem. B 101, 7821 (1997). doi
• A. Kovalenko and F. Hirata, "Three-dimensional density profiles of water in contact with a solute of arbitrary shape: a RISM approach," Chem. Phys. Lett. 290, 237 (1998). 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
• A. Kovalenko and F. Hirata, "Potentials of mean force of simple ions in ambient aqueous solution. I. Three-dimensional reference interaction site model approach," J. Chem. Phys. 112, 10391 (2000). 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. 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
• 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