ESM-RISM method¶

The ESM-RISM method is a hybrid solvation method that treats the electrode with density functional theory (DFT) and the ESM method, and the electrolyte with the integral equations of the reference interaction site model (RISM). It rewrites the RISM equation in the Laue representation, in which only the surface-normal direction is kept in real space, and connects the electrolyte to the bulk solution far from the electrode. Because both the electrode electrons and the solvent molecules and ions form grand-canonical ensembles, a single calculation describes a biased electrode together with the electric double layer that screens it. The method was introduced by Nishihara and Otani (2017), combined with the dielectrically consistent RISM by Hagiwara et al. (2022), and reformulated for efficiency by Hagiwara and Otani (2026).

This page uses atomic units (\(\hbar=e=m_\mathrm{e}=1\)) and \(\beta=1/k_\mathrm{B}T\). Solvent sites are labeled \(\alpha,\gamma\), with charge \(q_\gamma\) and bulk number density \(\rho_\gamma\).

Background and motivation¶

When a bias is applied to an electrode, an external potentiostat injects excess charge into it, and counter ions move from the bulk solution to the interface to screen that charge, forming an electric double layer. Describing this process requires a system in which both the number of electrode electrons and the number of ions in solution can change, that is, a grand-canonical ensemble.

First-principles molecular dynamics with an explicit electrolyte fixes the number of atoms, so it cannot add or remove counter ions flexibly, and converging thermodynamic quantities is very expensive. The 3D-RISM-SCF method, which combines 3D-RISM with DFT, obtains the solvent distribution efficiently from integral equations, but it treats electrostatics under three-dimensional periodic boundary conditions. Once the electrode carries a net charge, the periodic electrostatic potential diverges unless a uniform background charge is introduced. A periodically repeated cell also cannot represent the half-cell boundary condition, in which the electrolyte joins a semi-infinite bulk solution far from the electrode.

ESM-RISM removes both problems by treating surface-normal electrostatics with the open boundary conditions of the ESM method and by solving the RISM equation in real space along the same direction.

System setup and boundary conditions¶

Simulation cell¶

The \(z\) axis is normal to the surface, and periodic boundary conditions apply in the plane. In the typical half-cell model (vacuum/slab/solution, VSS), vacuum lies to the left of the electrode slab and the electrolyte to the right. The slab and the nearby solution sit inside the DFT unit cell; the solution continues through an expanded RISM region attached to the right edge of the DFT cell and joins the bulk solution far from the electrode. For molecules and ions in solution, the solution/ion/solution boundary condition places solution on both sides.

The electrons and ions treated by DFT (the solute) are the electrode and the adsorbates; the solvent sites treated by RISM belong to water molecules and ions. The number of solute electrons \(N\) can deviate from the neutral value \(N_0\) by \(\Delta N=N-N_0\), which gives the electrode an excess charge \(-\Delta N\).

Electrostatic potential from ESM¶

In the mixed representation of in-plane reciprocal vectors \(\boldsymbol{g}_{\parallel}\) and \(z\) (Laue representation), the Poisson equation reads

\[\begin{equation} \label{eq:poisson-laue} \partial_z\left[\epsilon(z)\,\partial_z V(\boldsymbol{g}_{\parallel},z)\right]-\epsilon(z)\,g_{\parallel}^2\,V(\boldsymbol{g}_{\parallel},z)=-4\pi\rho(\boldsymbol{g}_{\parallel},z). \end{equation}\]

ESM-RISM uses the vacuum/slab/vacuum ESM boundary condition (open boundaries with \(\epsilon(z)=1\) on both sides). The dielectric response of the solution enters explicitly as the RISM charge distribution, not as a dielectric continuum. The Green's function is then

\[\begin{equation} \label{eq:green-bc1} G(\boldsymbol{g}_{\parallel},z,z')= \begin{cases} \dfrac{2\pi}{g_{\parallel}}\,\mathrm{e}^{-g_{\parallel}|z-z'|} & (\boldsymbol{g}_{\parallel}\neq\boldsymbol{0}),\\[1ex] -2\pi|z-z'| & (\boldsymbol{g}_{\parallel}=\boldsymbol{0}), \end{cases} \end{equation}\]

and the electrostatic potential that an electron feels from the solute, \(v_\mathrm{es}(\boldsymbol{r})\) (the sum of the Hartree and local pseudopotential terms), is computed with this \(G\). Because the \(\boldsymbol{g}_{\parallel}=\boldsymbol{0}\) component is proportional to \(|z-z'|\), the potential of a net-charged system keeps changing linearly with \(z\) far away. The section on charge neutrality uses this property.

Laue representation of the RISM equation¶

From 3D-RISM to the Laue representation¶

In 3D-RISM, the total correlation function \(h_\gamma(\boldsymbol{r})\) of solvent site \(\gamma\) around the solute is related to the direct correlation function \(c_\alpha(\boldsymbol{r})\) by a convolution with the bulk solvent susceptibility \(\chi_{\alpha\gamma}\). ESM-RISM writes this relation in the Laue representation,

\[\begin{equation} \label{eq:lrism} 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'|), \end{equation}\]

and solves it (the Laue-represented RISM, or LRISM, equation). The in-plane convolution becomes a product in reciprocal space; only the convolution along \(z\) remains in real space. No periodicity along \(z\) is assumed, so the cell can have vacuum on one side and bulk solution on the other.

Solvent susceptibility¶

The solvent susceptibility is computed beforehand by a 1D-RISM calculation of the homogeneous, isotropic bulk solution. As a function of the radial wave number \(g=|\boldsymbol{g}|\),

\[\begin{equation} \label{eq:chi} \chi_{\alpha\gamma}(g)=\omega_{\alpha\gamma}(g)+\rho_\alpha h^\mathrm{VV}_{\alpha\gamma}(g), \end{equation}\]

where \(h^\mathrm{VV}_{\alpha\gamma}\) is the site–site total correlation function of the bulk solution and \(\omega_{\alpha\gamma}\) is the intramolecular correlation function. When \(\alpha\) and \(\gamma\) belong to the same molecule at distance \(l_{\alpha\gamma}\),

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

and it vanishes when they belong to different molecules. The Laue-represented susceptibility follows from a one-dimensional inverse Fourier transform along \(g_z\), with \(g=\sqrt{g_{\parallel}^2+g_z^2}\):

\[\begin{equation} \label{eq:chi-laue} \chi_{\alpha\gamma}(\boldsymbol{g}_{\parallel},|z-z'|)=\frac{1}{\pi}\int_0^\infty\mathrm{d}g_z\,\chi_{\alpha\gamma}\!\left(\sqrt{g_{\parallel}^2+g_z^2}\right)\cos\left[g_z(z-z')\right]. \end{equation}\]

All information on the bulk solution (molecular geometry and site–site correlations) enters through \(\chi_{\alpha\gamma}\).

Closure¶

Equation \(\eqref{eq:lrism}\) contains two unknowns, \(h_\gamma\) and \(c_\gamma\), so it is solved together with a closure relation that is local in real space. The hypernetted-chain (HNC) closure is

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

where \(g_\gamma=h_\gamma+1\) is the distribution function of the solvent site and \(u_\gamma\) is the solute–solvent interaction potential. In practice, the Kovalenko–Hirata (KH) closure, which linearizes the exponential in high-density regions, is used:

\[\begin{equation} \label{eq:kh} 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=-\beta u_\gamma+h_\gamma-c_\gamma. \end{equation}\]

Equation \(\eqref{eq:lrism}\) is evaluated in reciprocal space (\(\boldsymbol{g}_{\parallel}\)) and the closure in real space (\(\boldsymbol{r}_{\parallel}\)); the iteration switches between them with in-plane two-dimensional FFTs.

Solute–solvent interaction potential¶

Lennard-Jones and electrostatic terms¶

The interaction between solute atom \(I\) (at \(\boldsymbol{R}_I\)) and solvent site \(\gamma\) is the sum of a 12-6 Lennard-Jones (LJ) potential and an electrostatic term:

\[\begin{equation} \label{eq:u} 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{es}(\boldsymbol{r}). \end{equation}\]

The LJ parameters are combined with the Lorentz–Berthelot rules \(\sigma_{I\gamma}=(\sigma_I+\sigma_\gamma)/2\) and \(\varepsilon_{I\gamma}=\sqrt{\varepsilon_I\varepsilon_\gamma}\). Since \(v_\mathrm{es}\) is the potential energy of an electron, a site of charge \(q_\gamma\) feels the electrostatic energy \(-q_\gamma v_\mathrm{es}\). Because \(v_\mathrm{es}\) is obtained with the Green's function \(\eqref{eq:green-bc1}\), the field of a charged electrode reaches the solution without any artificial cancellation from periodic images. The LJ potential is truncated at a finite cutoff radius; this cutoff controls the accuracy of the short-range part, and some electrolytes require a larger value.

Short- and long-range separation of the direct correlation function¶

For numerical stability, the direct correlation function is split into short- and long-range parts:

\[\begin{equation} \label{eq:c-split} 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}). \end{equation}\]

The long-range part is the asymptotic form of \(c_\gamma\), namely the electrostatic term smoothed with width \(\tau\):

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

In ESM-RISM, ESM determines the far-field behavior of the Laue-represented \(c^\mathrm{L}_\gamma(\boldsymbol{g}_{\parallel},z)\). Because \(c^\mathrm{L}_\gamma\) is fixed by the solute electron density before the iteration, only \(c^\mathrm{S}_\gamma\) is solved iteratively. With \(u^\mathrm{S}_\gamma=u_\gamma-u^\mathrm{L}_\gamma\), the terms \(-\beta u^\mathrm{L}_\gamma\) and \(c^\mathrm{L}_\gamma\) cancel inside the exponent of the closure \(\eqref{eq:hnc}\), giving

\[\begin{equation} \label{eq:hnc-short} g_\gamma(\boldsymbol{r})=\exp\left[-\beta u^\mathrm{S}_\gamma(\boldsymbol{r})+h_\gamma(\boldsymbol{r})-c^\mathrm{S}_\gamma(\boldsymbol{r})\right]. \end{equation}\]

Wall potential against penetration into the electrode¶

Depending on the electrode morphology and the force-field parameters, the solvent distribution can leak into interstitial regions of the electrode and give unphysical results. To prevent this, a wall potential is added to \(u_\gamma\) in the closure. It is obtained by integrating the LJ repulsion from fictitious wall atoms of density \(\rho_\mathrm{wall}\) that fill \(z<z_\mathrm{r}\) over the plane and along \(z\):

\[\begin{equation} \label{eq:wall} u^\mathrm{wall}_\gamma(z)= \begin{cases} \infty & (z\le z_\mathrm{r}),\\[1ex] \dfrac{8\pi\rho_\mathrm{wall}\sqrt{\varepsilon_\mathrm{wall}\varepsilon_\gamma}}{90\,(z-z_\mathrm{r})^{9}}\left(\dfrac{\sigma_\mathrm{wall}+\sigma_\gamma}{2}\right)^{12} & (z>z_\mathrm{r}). \end{cases} \end{equation}\]

Far-field boundary conditions and connection to the bulk solution¶

Boundary conditions on the total correlation function¶

The laterally averaged total correlation function (the \(\boldsymbol{g}_{\parallel}=\boldsymbol{0}\) component) equals \(-1\) on the electrode side, where no solvent can enter, and vanishes far from the electrode:

\[\begin{equation} \label{eq:bc} h_\gamma(\boldsymbol{0},z)\to \begin{cases} -1 & (\text{electrode side}),\\ 0 & (z\to\infty), \end{cases} \qquad c^\mathrm{S}_\gamma(\boldsymbol{g}_{\parallel},z)\to 0\quad(z\to\infty). \end{equation}\]

The condition \(h_\gamma\to 0\) means that every site density returns to its bulk value \(\rho_\gamma\) far away. Because the bulk correlations enter through \(\chi_{\alpha\gamma}\) in Eq. \(\eqref{eq:chi}\), the far-field solution is exactly the bulk solution obtained by 1D-RISM. Through this connection, the interfacial solution exchanges sites with an infinitely large bulk reservoir, and the numbers of solvent molecules and ions are outputs of the calculation rather than inputs. This is the grand-canonical condition on the solution side.

The in-plane varying components \(h_\gamma(\boldsymbol{g}_{\parallel}\neq\boldsymbol{0},z)\) are localized near the electrode and vanish far away.

Expanded RISM region¶

Beyond the right edge of the DFT cell, the solution is treated in an expanded RISM region where only the lateral average matters. In this region only \(c^\mathrm{S}_\gamma(\boldsymbol{0},z)\) and \(h_\gamma(\boldsymbol{0},z)\) are solved, and the \(\boldsymbol{g}_{\parallel}\neq\boldsymbol{0}\) components are set to zero. The long-range part \(c^\mathrm{L}_\gamma(\boldsymbol{g}_{\parallel},z)\) is defined over the whole DFT cell plus expanded region so that its convolution is evaluated accurately. If the in-plane variation of \(c^\mathrm{S}_\gamma\) or \(h_\gamma\) has not decayed by the edge of the DFT cell, the calculation breaks down, so the solution region inside the DFT cell must be thick enough.

Stoichiometric constraint¶

Under the grand-canonical condition, the number of particles of each site is obtained by integrating the total correlation function,

\[\begin{equation} \label{eq:Nalpha} N_\alpha=S\rho_\alpha\int\mathrm{d}z\left[h_\alpha(\boldsymbol{0},z)+1\right], \end{equation}\]

where \(S\) is the lateral area. During the iteration, numerical error can make \(N_\alpha\) differ between sites of the same molecule and slow convergence. A correction with the weight function \(w_\alpha(z)=\mathrm{erfc}[(z-z_\mathrm{t})/\sigma_\mathrm{w}]/2\) is therefore added to the uncorrected total correlation function \(h^0_\alpha\):

\[\begin{equation} \label{eq:stoich} h_\alpha(\boldsymbol{g}_{\parallel},z)=h^0_\alpha(\boldsymbol{g}_{\parallel},z)+\eta_\alpha w_\alpha(z),\qquad \eta_\alpha=\frac{\tilde{N}_M-N^0_\alpha}{S\rho_\alpha\int\mathrm{d}z\,w_\alpha(z)}. \end{equation}\]

Here \(N^0_\alpha\) is the particle number computed from \(h^0_\alpha\), \(\tilde{N}_M\) is the average particle number over the sites of molecule \(M\), and \(z_\mathrm{t}\) is the position where \(h^0_\alpha(\boldsymbol{0},z)\) reaches zero.

Charge neutrality and the electric double layer¶

The charge density of the solution is

\[\begin{equation} \label{eq:rho-solv} \rho_\mathrm{solv}(\boldsymbol{r})=\sum_\gamma q_\gamma\rho_\gamma g_\gamma(\boldsymbol{r}). \end{equation}\]

Because the bulk solution is neutral (\(\sum_\gamma q_\gamma\rho_\gamma=0\)), \(\rho_\mathrm{solv}\) is nonzero only near the interface.

Let \(\bar{\rho}_\mathrm{tot}(z)\) be the laterally averaged total charge of solute and solution. The \(\boldsymbol{g}_{\parallel}=\boldsymbol{0}\) component of Eq. \(\eqref{eq:green-bc1}\) gives the laterally averaged total electrostatic potential

\[\begin{equation} \label{eq:phi0} \bar{\phi}(z)=-2\pi\int\mathrm{d}z'\,|z-z'|\,\bar{\rho}_\mathrm{tot}(z'), \end{equation}\]

whose slope outside the charge distribution is proportional to the total charge. A solution that returns to the bulk far away (\(h_\gamma\to 0\)) cannot sustain an electric field inside the bulk. The potential therefore becomes constant far away, at the inner potential \(\Phi_\mathrm{S}\) of the solution, only if the whole system is neutral:

\[\begin{equation} \label{eq:neutral} \Delta Q_\mathrm{RISM}\equiv\int\mathrm{d}\boldsymbol{r}\,\rho_\mathrm{solv}(\boldsymbol{r})=\Delta N. \end{equation}\]

The excess electrons \(\Delta N\) injected into the electrode are screened by the excess counter-ion charge \(\Delta Q_\mathrm{RISM}\) that accumulates at the interface; this is the electric double layer. Both the electrode charge and the counter charge are localized near the interface, so no uniform background charge is needed. Even with zero net charge, the double layer carries a dipole, which produces a potential difference between the vacuum side and the bulk solution. This potential difference is directly tied to the electrode potential.

Free energy and coupling to the Kohn–Sham equation¶

Helmholtz free energy¶

The total energy of the system is the Helmholtz free energy, the DFT total energy \(E_\mathrm{DFT}\) plus the excess chemical potential \(\Delta\mu_\mathrm{solv}\) of the solution:

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

With the KH closure, the solvation free energy has the closed form

\[\begin{equation} \label{eq:dmu} \Delta\mu_\mathrm{solv}=\frac{1}{\beta}\sum_\gamma\rho_\gamma\int\mathrm{d}\boldsymbol{r}\left[\frac{1}{2}h_\gamma^2(\boldsymbol{r})\,\Theta\!\left(-h_\gamma(\boldsymbol{r})\right)-c_\gamma(\boldsymbol{r})-\frac{1}{2}h_\gamma(\boldsymbol{r})\,c_\gamma(\boldsymbol{r})\right], \end{equation}\]

where \(\Theta\) is the step function. Dropping the first term gives the Gaussian fluctuation (GF) expression. The long-range part of \(c_\gamma\) does not vanish far away, but its contributions from the bulk region cancel because \(\sum_\gamma q_\gamma\rho_\gamma=0\).

Kohn–Sham equation¶

Varying Eq. \(\eqref{eq:A}\) with respect to the electron density \(n\) gives the Kohn–Sham equation including the solution:

\[\begin{equation} \label{eq:ks} \left[-\frac{\nabla^2}{2}+v_\mathrm{es}(\boldsymbol{r})+v_\mathrm{solv}(\boldsymbol{r})+v_\mathrm{xc}(\boldsymbol{r})\right]\psi_i(\boldsymbol{r})=\epsilon_i\psi_i(\boldsymbol{r}). \end{equation}\]

The solvation potential is the electrostatic potential of the solution charge, evaluated with the ESM Green's function:

\[\begin{equation} \label{eq:vsolv} v_\mathrm{solv}(\boldsymbol{r})=\frac{\delta\Delta\mu_\mathrm{solv}}{\delta n(\boldsymbol{r})}=-\int\mathrm{d}\boldsymbol{r}'\,G(\boldsymbol{r},\boldsymbol{r}')\,\rho_\mathrm{solv}(\boldsymbol{r}'). \end{equation}\]

Electrons are negatively charged, so positive counter ions attract them. In practice, the Kohn–Sham and LRISM equations are solved alternately until self-consistency: the electrons act on the solution through \(v_\mathrm{es}\), and the solution acts on the electrons through \(v_\mathrm{solv}\).

Janak's theorem and the electron chemical potential¶

Because the free-energy functional \(\eqref{eq:A}\) and the effective potential in Eq. \(\eqref{eq:ks}\) are variationally consistent, Janak's theorem also holds in ESM-RISM:

\[\begin{equation} \label{eq:janak} \frac{\partial A}{\partial f_i}=\epsilon_i,\qquad \mu_\mathrm{e}=\frac{\partial A}{\partial N}, \end{equation}\]

where \(f_i\) are occupation numbers. The Fermi level measured from the inner potential \(\Phi_\mathrm{S}\) of the solution therefore equals the change in free energy per added electron.

Energy reference¶

Under the vacuum/slab/vacuum ESM boundary condition, the origin of the electrostatic potential is the midpoint of its values at the two cell edges. To measure the electrode potential from the bulk solution, \(v_\mathrm{es}\) is shifted by \(v_0\) so that the electrostatic potential of the bulk solution is zero. This shift changes both the electrostatic energy of a charged system and \(\Delta\mu_\mathrm{solv}\); to keep Janak's theorem, the free energy is corrected to

\[\begin{equation} \label{eq:A-shift} A=T_\mathrm{s}+E_\mathrm{es}+E_\mathrm{xc}+\Delta\mu_\mathrm{solv}-\frac{v_0\,\Delta N}{2}, \end{equation}\]

where \(E_\mathrm{es}\) is the sum of the Hartree and ionic electrostatic energies. The correction appears only for charged systems (\(\Delta N\neq 0\)).

Grand potential and the constant-μ method¶

To compare systems with different numbers of electrons, ESM-RISM uses the grand potential obtained by a Legendre transform with respect to the electrons,

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

The solution is already in equilibrium with the bulk through the boundary conditions of the LRISM equation, so \(\Omega\) is the thermodynamic potential of a system open to both electrons and solution particles.

Instead of scanning \(\Delta N\) and computing \(\Omega\) and \(\mu_\mathrm{e}\), one can fix \(\mu_\mathrm{e}\) at a target value with the constant-μ method and let the electron number adjust. In ESM-RISM the counter charge in the solution follows the excess electrode charge self-consistently, so combining it with the constant-μ method allows structure optimizations and reaction-path calculations at a fixed electrode potential.

Combination with the dielectrically consistent RISM¶

The conventional 1D-RISM (extended RISM, XRISM) systematically underestimates the relative permittivity of solutions of dipolar solvent molecules. The dielectrically consistent RISM (DRISM) fixes this by taking the experimental permittivity \(\varepsilon\) as a parameter and adding a correction function \(\xi_{\alpha\gamma}\). In reciprocal-space matrix form, the 1D-RISM equation becomes

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

where \(\rho\) is the diagonal matrix of site number densities. With each molecule oriented so that its dipole lies along \(z\), the correction function between polar-solvent sites is

\[\begin{equation} \label{eq:xi} \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), \end{equation}\]

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

Here \(j_0\) and \(j_1\) are spherical Bessel functions, \(d^{x,y,z}_\alpha\) are site positions measured from the charge center of the molecule, \(\rho_\mathrm{tot}\) is the total number density, \(\rho_\mathrm{pol}\) and \(\mu\) are the number density and dipole moment of the polar molecules, and \(a\) is a smoothing width (about 0.5 Å for aqueous solutions).

The DRISM solution gives the solvent susceptibility in the same form as Eq. \(\eqref{eq:chi}\):

\[\begin{equation} \label{eq:chi-drism} \chi=\omega'+\rho h'=\omega+\rho h. \end{equation}\]

ESM-RISM receives information on the solution only through \(\chi_{\alpha\gamma}\), so the DRISM improvement carries over to interface calculations without changing the LRISM equation. For the Pt(111)/1 M HCl(aq) interface, DRISM brings the potential of zero charge and the double-layer capacitance closer to experiment (Hagiwara et al. 2022).

Numerical solution of the LRISM equation¶

Treatment of the intramolecular correlation function¶

Inverse-transforming \(\omega_{\alpha\gamma}\) directly in Eq. \(\eqref{eq:chi-laue}\) produces numerical noise, because in real space \(\omega_{\alpha\gamma}(r)=\delta(r-l_{\alpha\gamma})/(4\pi r^2)\) is a delta function. The original paper replaced the delta function with a Gaussian of width \(\kappa\):

\[\begin{equation} \label{eq:omega-gauss} \omega_{\alpha\gamma}(g)=\frac{\sin(g\,l_{\alpha\gamma})}{g\,l_{\alpha\gamma}}\exp\left(-\frac{g^2\kappa^2}{4}\right). \end{equation}\]

This approximation gives the 1D-RISM intramolecular correlation an artificial width, so ESM-RISM results deviate from those of 3D-RISM-SCF. Hagiwara et al. (2022) also introduced a damping scheme that leaves \(\omega_{\alpha\gamma}\) unchanged, truncates the integral at \(g_\mathrm{cut}\), and multiplies it by a damping function:

\[\begin{equation} \label{eq:damping} \chi^{\omega}_{\alpha\gamma}(\boldsymbol{g}_{\parallel},|z-z'|)=\frac{1}{\pi}\int_0^{g_\mathrm{cut}}\mathrm{d}g_z\,\omega_{\alpha\gamma}(g)\,W(g_z)\cos\left[g_z(z-z')\right],\qquad W(g_z)=1-3\left(\frac{g_z}{g_\mathrm{cut}}\right)^2+2\left(\frac{g_z}{g_\mathrm{cut}}\right)^3, \end{equation}\]

where \(\chi^{\omega}_{\alpha\gamma}\) denotes the \(\omega_{\alpha\gamma}\) contribution to Eq. \(\eqref{eq:chi-laue}\).

Separating the in-plane and laterally averaged components¶

Solving Eq. \(\eqref{eq:lrism}\) as written requires a real-space convolution along \(z\) for every \(\boldsymbol{g}_{\parallel}\), which costs \(O(N^{4/3})\) for \(N\) total grid points. Hagiwara and Otani (2026) removed this bottleneck by splitting the LRISM equation into its \(\boldsymbol{g}_{\parallel}\neq\boldsymbol{0}\) and \(\boldsymbol{g}_{\parallel}=\boldsymbol{0}\) components.

The \(\boldsymbol{g}_{\parallel}\neq\boldsymbol{0}\) components are localized near the electrode, so they can also be Fourier transformed along \(z\) and convolved as products in reciprocal space, as in 3D-RISM. The cost is that of a three-dimensional FFT, \(O(N\log N)\).

The \(\boldsymbol{g}_{\parallel}=\boldsymbol{0}\) component changes from \(-1\) to 0, so it is solved in real space without Fourier transforms. Substituting Eq. \(\eqref{eq:c-split}\) gives

\[\begin{equation} \label{eq:lrism-g0} h_\gamma(\boldsymbol{0},z)=\sum_\alpha\int\mathrm{d}z'\,c^\mathrm{S}_\alpha(\boldsymbol{0},z')\,\chi_{\alpha\gamma}(\boldsymbol{0},|z-z'|)+\sum_\alpha\int\mathrm{d}z'\,c^\mathrm{L}_\alpha(\boldsymbol{0},z')\,\chi_{\alpha\gamma}(\boldsymbol{0},|z-z'|). \end{equation}\]

For the second (long-range) term, the electrostatic potential is split into the periodic solution and the ESM correction:

\[\begin{equation} \label{eq:v-split} v_\mathrm{es}(\boldsymbol{r})=v^\mathrm{PBC}_\mathrm{es}(\boldsymbol{r})+v^\mathrm{ESM}_\mathrm{es}(\boldsymbol{r}). \end{equation}\]

The periodic part \(v^\mathrm{PBC}_\mathrm{es}\) is convolved in reciprocal space. Because \(v^\mathrm{ESM}_\mathrm{es}(\boldsymbol{0},z)\) is quadratic in \(z\), its contribution is obtained analytically by integration in spherical coordinates. The long-range term needs to be computed only once, before the iteration.

For the first (short-range) term, \(\chi_{\alpha\gamma}\) is split into \(\omega_{\alpha\gamma}\) and \(\rho_\alpha h^\mathrm{VV}_{\alpha\gamma}\). The in-plane integral of the intramolecular correlation function is analytic:

\[\begin{equation} \label{eq:omega-g0} \omega_{\alpha\gamma}(\boldsymbol{0},|z-z'|)= \begin{cases} \dfrac{1}{2l_{\alpha\gamma}} & (|z-z'|<l_{\alpha\gamma}),\\[1ex] \dfrac{1}{4l_{\alpha\gamma}} & (|z-z'|=l_{\alpha\gamma}),\\[1ex] 0 & (|z-z'|>l_{\alpha\gamma}), \end{cases} \end{equation}\]

so neither the Gaussian nor the damping approximation is needed (for \(l_{\alpha\gamma}=0\) it reduces to \(\delta(z-z')\)). The in-plane integral of \(h^\mathrm{VV}_{\alpha\gamma}\) is tabulated beforehand in cylindrical coordinates. This real-space convolution costs \(O(N^{2/3})\).

The dominant cost of the iteration thus becomes the \(O(N\log N)\) of the \(\boldsymbol{g}_{\parallel}\neq\boldsymbol{0}\) components. The hydration free energy of a water molecule in pure water then agrees with 3D-RISM-SCF to within \(1.6\times10^{-4}\) kcal/mol (compared with \(+0.66\) kcal/mol for the Gaussian approximation and \(-0.15\) kcal/mol with the window function).

Physical quantities¶

Electrode potential¶

With the electron chemical potential \(\mu_\mathrm{e}\) measured from the inner potential \(\Phi_\mathrm{S}\) of the solution, the electrode potential is

\[\begin{equation} \label{eq:U} U[\text{vs.}\ \Phi_\mathrm{S}]=-\mu_\mathrm{e}/e. \end{equation}\]

Its value at zero excess charge is the potential of zero charge (PZC), \(U_\mathrm{PZC}\). The standard hydrogen electrode (SHE) potential is obtained from the free-energy change \(\Delta G_\mathrm{SHE}\) of the reaction \(\mathrm{H_2O(aq)}+\tfrac{1}{2}\mathrm{H_2(g)}\rightleftharpoons\mathrm{H_3O^+(aq)}+\mathrm{e^-(M)}\) as \(U_\mathrm{SHE}=-\Delta G_\mathrm{SHE}/nF\), computed relative to \(\Phi_\mathrm{S}\) by ESM-RISM with the same solution parameters. Then

\[\begin{equation} \label{eq:U-SHE} U[\text{vs. SHE}]=U[\text{vs.}\ \Phi_\mathrm{S}]-U_\mathrm{SHE}[\text{vs.}\ \Phi_\mathrm{S}] \end{equation}\]

gives the electrode potential on the same scale as experiment.

Double-layer capacitance¶

Viewed as a function of the electrode potential, the grand potential \(\eqref{eq:omega-grand}\) satisfies, by \(\partial A/\partial N=\mu_\mathrm{e}\),

\[\begin{equation} \label{eq:lippmann} \frac{\partial\Omega}{\partial U}=e\,\Delta N,\qquad C_\mathrm{dl}=-\frac{\partial^2\Omega}{\partial U^2}. \end{equation}\]

The first relation is the Lippmann equation for the electrode charge \(-e\Delta N\). \(\Omega\) is close to a downward parabola with its maximum at the PZC, and its curvature gives the double-layer capacitance \(C_\mathrm{dl}\).

Applications¶

For the Pt(111)/1 M HCl(aq) interface, the PZC, the SHE potential, and the double-layer capacitance near the PZC agree with experiment (Hagiwara et al. 2022). With the reformulated method, the PZC is \(+0.23\) V vs. SHE, and the computed SHE potential converted to the absolute scale is 4.06–4.37 V, consistent with the IUPAC recommended value of 4.44 V (Hagiwara and Otani 2026). The charge and electrostatic-potential profiles of the solution and the potential dependence of the grand potential of adsorbed species also give access to potential-dependent adsorption sites and equilibrium potentials of reactions.

References¶

M. Otani and O. Sugino, "First-principles calculations of charged surfaces and interfaces: A plane-wave nonrepeated slab approach," Phys. Rev. B 73, 115407 (2006). 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. Materials 6, 093802 (2022). doi
S. Hagiwara and M. Otani, "An Efficient Implementation of the ESM-RISM Method for Simulating Electrode/Electrolyte Interfaces," J. Chem. Theory Comput. 22, 3633 (2026). 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
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
J. Haruyama, T. Ikeshoji, and M. Otani, "Electrode potential from density functional theory calculations combined with implicit solvation theory," Phys. Rev. Materials 2, 095801 (2018). doi
N. Bonnet, T. Morishita, O. Sugino, and M. Otani, "First-Principles Molecular Dynamics at a Constant Electrode Potential," Phys. Rev. Lett. 109, 266101 (2012). doi