ESM-RISM法

ESM-RISM法は、電極を密度汎関数理論（DFT）とESM法で、電解液を参照相互作用サイト模型（RISM）の積分方程式で扱うハイブリッド溶媒和法です。RISM方程式を表面垂直方向だけ実空間で扱うLaue表示に書き換え、電解液を電極から遠方のバルク溶液に接続します。電極の電子と溶液中の分子・イオンをともにグランドカノニカル集団として扱えるため、バイアス電位を印加した電極と、それを遮蔽する電気二重層を一つの計算で記述できます。原論文はNishihara and Otani (2017)、誘電的に整合なRISMとの結合はHagiwara et al. (2022)、Laue表示RISM方程式の再定式化と効率化はHagiwara and Otani (2026)です。

本ページでは原子単位（\(\hbar=e=m_\mathrm{e}=1\)）を用い、\(\beta=1/k_\mathrm{B}T\) とします。溶媒サイトを \(\alpha,\gamma\)、その電荷を \(q_\gamma\)、バルク数密度を \(\rho_\gamma\) と書きます。

背景と動機

電極電位を印加すると、外部のポテンショスタットから電極に余剰電荷が注入され、それを遮蔽するカウンターイオンが溶液バルクから界面に集まって電気二重層が形成されます。この過程を扱うには、電極の電子数と溶液中のイオン数がともに変化する系、すなわちグランドカノニカル集団を扱う必要があります。

第一原理分子動力学法で電解液を陽に扱う方法は、原子数が固定されるためカウンターイオンを柔軟に出し入れできず、熱力学量の収束にも膨大な計算量を要します。3D-RISM法をDFTと組み合わせた3D-RISM-SCF法は溶液分布を積分方程式で効率よく求めますが、静電相互作用を三次元周期境界条件で扱います。電極が帯電すると、周期境界条件の下では静電ポテンシャルが発散するため、一様な背景電荷を導入しない限り計算できません。また、周期的に繰り返されたセルでは、電極から十分遠方で半無限のバルク溶液に接続する、という半電池の境界条件を表現できません。

ESM-RISM法は、表面垂直方向の静電相互作用をESM法の開放境界条件で扱い、RISM方程式も同じ方向について実空間で解くことで、これらの問題を解決します。

系の設定と境界条件

計算セル

表面垂直方向を \(z\) 軸にとり、面内方向には周期境界条件を課します。典型的な半電池模型（vacuum/slab/solution、VSS）では、電極スラブの左側を真空、右側を電解液とします。電極スラブとその近傍の溶液はDFTのユニットセル内に置き、溶液はDFTセルの右端からさらに拡張したRISM領域を経て、電極から十分遠方でバルク溶液に接続します。溶液中の分子やイオンを扱うときは、両側を溶液とする境界条件（solution/ion/solution）も用います。

DFTで扱う電子とイオン（以下、溶質）は電極と吸着種、RISMで扱う溶媒サイトは水分子やイオンを構成するサイトです。溶質の電子数 \(N\) は電気的に中性な系の電子数 \(N_0\) から \(\Delta N=N-N_0\) だけずらせるので、電極に余剰電荷 \(-\Delta N\) を与えられます。

ESMによる静電ポテンシャル

面内方向の逆格子ベクトル \(\boldsymbol{g}_{\parallel}\) と \(z\) の混合表示（Laue表示）で、ポアソン方程式は

\[\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法では、ESMの境界条件として真空/スラブ/真空（両側とも \(\epsilon(z)=1\) の開放境界条件）を用います。溶液の誘電応答は連続誘電体ではなく、RISMで求めた溶液の電荷分布として陽に入るためです。このときグリーン関数は

\[\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}\]

で、電子が溶質から感じる静電ポテンシャル \(v_\mathrm{es}(\boldsymbol{r})\)（Hartree項と局所擬ポテンシャル項の和）はこの \(G\) を用いて求めます。\(\boldsymbol{g}_{\parallel}=\boldsymbol{0}\) 成分が \(|z-z'|\) に比例するため、系全体が帯電していると静電ポテンシャルは遠方で \(z\) に比例して変化し続けます。この性質は後述の電荷中性の議論で使います。

Laue表示のRISM方程式

3D-RISMからLaue表示へ

3D-RISM法では、溶質周りの溶媒サイト \(\gamma\) の全相関関数 \(h_\gamma(\boldsymbol{r})\) と直接相関関数 \(c_\alpha(\boldsymbol{r})\) を、バルク溶液の溶媒感受率 \(\chi_{\alpha\gamma}\) を介した畳み込みで結びます。ESM-RISM法ではこれをLaue表示で書き、

\[\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}\]

を解きます（Laue表示RISM、LRISM方程式）。面内方向は逆空間で積となり、\(z\) 方向だけが実空間の畳み込みとして残ります。このため、\(z\) 方向に周期性を仮定せずに、片側が真空でもう片側がバルク溶液という非対称な境界条件を課せます。

溶媒感受率

溶媒感受率は、一様等方なバルク溶液に対する1D-RISM計算からあらかじめ求めます。逆空間の動径 \(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}\]

と書けます。\(h^\mathrm{VV}_{\alpha\gamma}\) はバルク溶液のサイト間全相関関数、\(\omega_{\alpha\gamma}\) は分子内相関関数です。\(\alpha\) と \(\gamma\) が同じ分子に属し、その距離が \(l_{\alpha\gamma}\) のとき

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

で、異なる分子に属するときは0です。Laue表示の感受率は、\(g=\sqrt{g_{\parallel}^2+g_z^2}\) として \(g_z\) 方向に一次元逆フーリエ変換して得ます。

\[\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}\]

バルク溶液の構造（分子形状とサイト間相関）はすべて \(\chi_{\alpha\gamma}\) に含まれます。

クロージャ

式\(\eqref{eq:lrism}\)は未知関数 \(h_\gamma\) と \(c_\gamma\) を含むので、実空間で局所的なクロージャ方程式と連立させて解きます。hypernetted-chain（HNC）クロージャは

\[\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}\]

で、\(g_\gamma=h_\gamma+1\) は溶媒サイトの分布関数、\(u_\gamma\) は溶質と溶媒サイトの相互作用ポテンシャルです。実際の計算では、HNCの指数関数を密度の高い領域で線形化したKovalenko–Hirata（KH）クロージャ

\[\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}\]

を用います。式\(\eqref{eq:lrism}\)は逆空間（\(\boldsymbol{g}_{\parallel}\)）、クロージャは実空間（\(\boldsymbol{r}_{\parallel}\)）で評価し、面内の二次元高速フーリエ変換で行き来しながら反復します。

溶質–溶媒相互作用ポテンシャル

Lennard-Jones項と静電項

溶質原子 \(I\)（位置 \(\boldsymbol{R}_I\)）と溶媒サイト \(\gamma\) の相互作用は、12-6型Lennard-Jones（LJ）ポテンシャルと静電項の和で表します。

\[\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}\]

LJパラメータはLorentz–Berthelot則 \(\sigma_{I\gamma}=(\sigma_I+\sigma_\gamma)/2\)、\(\varepsilon_{I\gamma}=\sqrt{\varepsilon_I\varepsilon_\gamma}\) で組み合わせます。\(v_\mathrm{es}\) は電子に対するポテンシャルエネルギーなので、電荷 \(q_\gamma\) のサイトが感じる静電エネルギーは \(-q_\gamma v_\mathrm{es}\) です。\(v_\mathrm{es}\) を式\(\eqref{eq:green-bc1}\)のグリーン関数で求めるので、帯電した電極が作る静電場は、周期境界条件による人工的な打ち消しなしに溶液へ伝わります。LJポテンシャルは有限のカットオフ半径で打ち切ります。このカットオフは短距離成分の精度を左右するので、電解液によっては大きめにとる必要があります。

直接相関関数の短距離・長距離分割

数値的な安定性のため、直接相関関数を短距離成分と長距離成分に分けます。

\[\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}\]

長距離成分は \(c_\gamma\) の漸近形で、静電項を幅 \(\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}\]

ESM-RISM法では、このLaue表示 \(c^\mathrm{L}_\gamma(\boldsymbol{g}_{\parallel},z)\) の遠方での振る舞いをESMで決めます。\(c^\mathrm{L}_\gamma\) は溶質の電子密度から反復前に一度だけ決まるので、反復で求めるのは \(c^\mathrm{S}_\gamma\) だけです。\(u^\mathrm{S}_\gamma=u_\gamma-u^\mathrm{L}_\gamma\) とすると、クロージャ式\(\eqref{eq:hnc}\)の指数の中では \(-\beta u^\mathrm{L}_\gamma\) と \(c^\mathrm{L}_\gamma\) が打ち消し合い、

\[\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}\]

となります。

電極内部への浸透を防ぐ壁ポテンシャル

電極の形状や力場パラメータによっては、溶液の分布関数が電極内部の隙間に入り込み、非物理的な結果を与えます。これを防ぐため、\(z<z_\mathrm{r}\) を占める密度 \(\rho_\mathrm{wall}\) の仮想的な壁原子とのLJ斥力項を面内と \(z\) 方向に積分した壁ポテンシャルをクロージャの \(u_\gamma\) に加えます。

\[\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}\]

遠方の境界条件とバルク溶液への接続

全相関関数の境界条件

全相関関数の面内平均（\(\boldsymbol{g}_{\parallel}=\boldsymbol{0}\) 成分）は、溶媒が入り込めない電極側で \(-1\)、電極から十分遠方で0になります。

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

\(h_\gamma\to 0\) は、遠方で各サイトの密度がバルク密度 \(\rho_\gamma\) に戻ることを意味します。バルクの相関は式\(\eqref{eq:chi}\)の \(\chi_{\alpha\gamma}\) を通じて入っているので、遠方の溶液は1D-RISMで求めたバルク溶液そのものになります。この接続によって、界面付近の溶液は無限に大きなバルク溶液とサイトの交換ができ、溶媒分子やイオンの数は入力ではなく計算結果として決まります。これが溶液側のグランドカノニカル条件です。

一方、面内で変化する成分 \(h_\gamma(\boldsymbol{g}_{\parallel}\neq\boldsymbol{0},z)\) は電極近傍にだけ局在し、遠方では0になります。

拡張RISM領域

DFTセルの右端から先は、溶液の面内平均だけが意味を持つ拡張RISM領域として扱います。この領域では \(c^\mathrm{S}_\gamma(\boldsymbol{0},z)\) と \(h_\gamma(\boldsymbol{0},z)\) だけを解き、\(\boldsymbol{g}_{\parallel}\neq\boldsymbol{0}\) 成分を0とします。長距離成分 \(c^\mathrm{L}_\gamma(\boldsymbol{g}_{\parallel},z)\) は、畳み込みを精度よく評価するため、DFTセルと拡張領域の全体で定義します。\(c^\mathrm{S}_\gamma\) や \(h_\gamma\) の面内変化がDFTセルの端までに十分減衰していないと計算が破綻するので、DFTセル内の溶液領域は十分厚くとる必要があります。

化学量論の拘束

グランドカノニカル条件では、各サイトの粒子数を全相関関数の積分から

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

と求めます（\(S\) は面内の断面積）。反復の途中では、同一分子の異なるサイトの \(N_\alpha\) が数値誤差でそろわず、収束を妨げることがあります。そこで、補正前の全相関関数 \(h^0_\alpha\) に重み関数 \(w_\alpha(z)=\mathrm{erfc}[(z-z_\mathrm{t})/\sigma_\mathrm{w}]/2\) を掛けた補正を加え、

\[\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}\]

とします。\(N^0_\alpha\) は \(h^0_\alpha\) から求めた粒子数、\(\tilde{N}_M\) は分子 \(M\) を構成するサイトの粒子数の平均、\(z_\mathrm{t}\) は \(h^0_\alpha(\boldsymbol{0},z)\) が0になる位置です。

電荷中性と電気二重層

溶液の電荷密度は

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

です。バルク溶液は中性（\(\sum_\gamma q_\gamma\rho_\gamma=0\)）なので、\(\rho_\mathrm{solv}\) は界面近傍にだけ現れます。

溶質と溶液を合わせた全電荷の面内平均を \(\bar{\rho}_\mathrm{tot}(z)\) とすると、式\(\eqref{eq:green-bc1}\)の \(\boldsymbol{g}_{\parallel}=\boldsymbol{0}\) 成分から、全静電ポテンシャルの面内平均は

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

で、電荷分布の外側での傾きは全電荷に比例します。溶液が遠方でバルクに戻る（\(h_\gamma\to 0\)）ことと、バルク溶液中に電場が残ることは両立しません。したがって遠方で静電ポテンシャルが一定値（溶液の内部ポテンシャル \(\Phi_\mathrm{S}\)）になるには、系全体が中性である必要があります。

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

電極に注入した余剰電子 \(\Delta N\) は、界面に集まる過剰なカウンターイオンの電荷 \(\Delta Q_\mathrm{RISM}\) で遮蔽されます。これが電気二重層です。電極の電荷とカウンターチャージがともに界面近傍に局在するので、一様な背景電荷を導入する必要はありません。全電荷が0でも電気二重層は双極子を持つので、真空側と溶液バルクの間には電位差が生じます。この電位差が電極電位と直接結びつきます。

自由エネルギーとKohn–Sham方程式への結合

ヘルムホルツ自由エネルギー

系の全エネルギーは、DFTの全エネルギー \(E_\mathrm{DFT}\) に溶液による過剰化学ポテンシャル \(\Delta\mu_\mathrm{solv}\) を加えたヘルムホルツ自由エネルギーです。

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

KHクロージャを用いると、溶媒和自由エネルギーは解析的に

\[\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}\]

と書けます（\(\Theta\) は階段関数）。第1項を除いたものがガウスゆらぎ（GF）近似の表式です。\(c_\gamma\) の長距離成分は遠方でも0になりませんが、\(\sum_\gamma q_\gamma\rho_\gamma=0\) なのでバルク領域の寄与は打ち消し合います。

Kohn–Sham方程式

式\(\eqref{eq:A}\)を電子密度 \(n\) で変分すると、溶液の効果を含むKohn–Sham方程式

\[\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}\]

が得られます。溶媒和ポテンシャルは溶液の電荷分布がつくる静電ポテンシャルで、ESMのグリーン関数を用いて

\[\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}\]

と書けます。電子の電荷が負なので、正のカウンターイオンは電子を引きつけます。実際の計算では、Kohn–Sham方程式とLRISM方程式を交互に解き、\(v_\mathrm{es}\) を通じて電子から溶液へ、\(v_\mathrm{solv}\) を通じて溶液から電子へ相互作用を伝えて自己無撞着解を求めます。

Janakの定理と電子の化学ポテンシャル

自由エネルギー汎関数\(\eqref{eq:A}\)と有効ポテンシャル\(\eqref{eq:ks}\)が変分的に整合しているので、ESM-RISM法でもJanakの定理

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

が成り立ちます（\(f_i\) は占有数）。つまり、Fermi準位を溶液の内部ポテンシャル \(\Phi_\mathrm{S}\) から測った値が、電子数を変えたときの自由エネルギー変化に一致します。

エネルギーの基準

ESMの真空/スラブ/真空境界条件では、静電ポテンシャルの原点はセル両端の値の中点にとられます。電極電位を溶液バルク基準で議論するため、溶液バルクの静電ポテンシャルが0になるように \(v_\mathrm{es}\) を \(v_0\) だけずらします。このずらしは帯電した系の静電エネルギーと \(\Delta\mu_\mathrm{solv}\) の両方を変えるので、Janakの定理を保つには自由エネルギーを

\[\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}\]

と補正します。\(E_\mathrm{es}\) はHartreeエネルギーとイオンの静電エネルギーの和です。補正項は帯電した系（\(\Delta N\neq 0\)）でだけ現れます。

グランドポテンシャルとConstant-μ法

電子数の異なる系を比較するには、電子についてルジャンドル変換したグランドポテンシャル

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

を用います。溶液側はLRISM方程式の境界条件によってすでにバルク溶液と平衡にあるので、\(\Omega\) は電子と溶液粒子の両方について開いた系の熱力学ポテンシャルです。

\(\Delta N\) を変えて \(\Omega\) と \(\mu_\mathrm{e}\) を求める代わりに、Constant-μ法で \(\mu_\mathrm{e}\) を目標値に固定し、電子数を自動的に調整することもできます。ESM-RISM法では電極の余剰電荷に応じて溶液のカウンターチャージが自己無撞着に決まるので、Constant-μ法と組み合わせると、電極電位を一定に保ったまま構造最適化や反応経路の計算を行えます。

誘電的に整合なRISMとの結合

通常の1D-RISM（拡張RISM、XRISM）は、双極子を持つ溶媒分子からなる溶液の比誘電率を系統的に過小評価します。誘電的に整合なRISM（dielectrically consistent RISM、DRISM）は、実験の比誘電率 \(\varepsilon\) をパラメータとして与え、補正関数 \(\xi_{\alpha\gamma}\) を加えてこれを改善します。1D-RISM方程式を逆空間の行列形式で

\[\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}\]

と書きます（\(\rho\) はサイト数密度の対角行列）。分子が \(z\) 軸方向に双極子を持つように置いたとき、極性溶媒のサイト間の補正関数は

\[\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}\]

です。\(j_0,j_1\) は球ベッセル関数、\(d^{x,y,z}_\alpha\) は分子の電荷中心から測ったサイト位置、\(\rho_\mathrm{tot}\) は全数密度、\(\rho_\mathrm{pol}\) と \(\mu\) は極性分子の数密度と双極子モーメント、\(a\) は平滑化幅（水溶液では0.5 Å程度）です。

DRISMの解からも溶媒感受率は

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

で、式\(\eqref{eq:chi}\)と同じ形で得られます。ESM-RISM法は溶液の情報を \(\chi_{\alpha\gamma}\) だけを通じて受け取るので、LRISM方程式に手を加えずにDRISMの改善を界面計算に取り込めます。Pt(111)/1 M HCl水溶液界面では、DRISMを用いるとゼロ電荷電位と電気二重層容量が実験値に近づきます（Hagiwara et al. 2022）。

LRISM方程式の数値解法

分子内相関関数の扱い

式\(\eqref{eq:chi-laue}\)で \(\omega_{\alpha\gamma}\) をそのまま逆フーリエ変換すると、実空間の \(\omega_{\alpha\gamma}(r)=\delta(r-l_{\alpha\gamma})/(4\pi r^2)\) がデルタ関数なので数値ノイズが出ます。原論文では、デルタ関数を幅 \(\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}\]

この近似は1D-RISMの分子内相関関数に人工的な幅を与えるので、ESM-RISM法の結果が3D-RISM-SCF法の結果とずれます。Hagiwara et al. (2022)は、\(\omega_{\alpha\gamma}\) を変えずに、積分範囲を \(g_\mathrm{cut}\) までに制限して減衰関数を掛ける方法（damping scheme）も導入しました。

\[\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}\]

ここで \(\chi^{\omega}_{\alpha\gamma}\) は式\(\eqref{eq:chi-laue}\)のうち \(\omega_{\alpha\gamma}\) の寄与です。

面内成分と面内平均成分の分離

式\(\eqref{eq:lrism}\)をそのまま解くと、各 \(\boldsymbol{g}_{\parallel}\) について \(z\) 方向の畳み込みを実空間で行うので、計算量は \(O(N^{4/3})\)（\(N\) は全格子点数）になります。Hagiwara and Otani (2026)は、LRISM方程式を \(\boldsymbol{g}_{\parallel}\neq\boldsymbol{0}\) 成分と \(\boldsymbol{g}_{\parallel}=\boldsymbol{0}\) 成分に分けてこれを解消しました。

\(\boldsymbol{g}_{\parallel}\neq\boldsymbol{0}\) 成分は電極近傍に局在するので、\(z\) 方向にもフーリエ変換でき、3D-RISMと同じく逆空間の積で畳み込めます。計算量は三次元高速フーリエ変換の \(O(N\log N)\) です。

\(\boldsymbol{g}_{\parallel}=\boldsymbol{0}\) 成分は \(-1\) から0へ変化するので、フーリエ変換を使わずに実空間で解きます。式\(\eqref{eq:c-split}\)を代入すると

\[\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}\]

となります。第2項（長距離項）では、静電ポテンシャルを周期境界条件の解とESMの補正に分けます。

\[\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}\]

周期的な \(v^\mathrm{PBC}_\mathrm{es}\) の寄与は逆空間で畳み込めます。\(v^\mathrm{ESM}_\mathrm{es}(\boldsymbol{0},z)\) は \(z\) の二次関数なので、その寄与は球座標の積分で解析的に求まります。長距離項は反復の前に一度計算すれば十分です。

第1項（短距離項）では、\(\chi_{\alpha\gamma}\) を \(\omega_{\alpha\gamma}\) と \(\rho_\alpha h^\mathrm{VV}_{\alpha\gamma}\) に分けます。分子内相関関数の面内積分は解析的に

\[\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}\]

となるので、ガウス関数や減衰関数による近似は不要です（\(l_{\alpha\gamma}=0\) では \(\delta(z-z')\)）。\(h^\mathrm{VV}_{\alpha\gamma}\) の面内積分は円柱座標であらかじめ用意できます。この実空間の畳み込みの計算量は \(O(N^{2/3})\) です。

この結果、反復の主な計算量は \(\boldsymbol{g}_{\parallel}\neq\boldsymbol{0}\) 成分の \(O(N\log N)\) になります。純水中の水分子の水和自由エネルギーは、3D-RISM-SCF法との差が \(1.6\times10^{-4}\) kcal/molまで小さくなります（ガウス近似では \(+0.66\) kcal/mol、Window関数では \(-0.15\) kcal/mol）。

得られる物理量

電極電位

電子の化学ポテンシャル \(\mu_\mathrm{e}\) を溶液の内部ポテンシャル \(\Phi_\mathrm{S}\) から測ると、電極電位は

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

で与えられます。余剰電荷0での値がゼロ電荷電位（PZC）\(U_\mathrm{PZC}\) です。標準水素電極（SHE）の電位は、反応 \(\mathrm{H_2O(aq)}+\tfrac{1}{2}\mathrm{H_2(g)}\rightleftharpoons\mathrm{H_3O^+(aq)}+\mathrm{e^-(M)}\) の自由エネルギー変化 \(\Delta G_\mathrm{SHE}\) から \(U_\mathrm{SHE}=-\Delta G_\mathrm{SHE}/nF\) として、同じ溶液パラメータのESM-RISM計算で \(\Phi_\mathrm{S}\) 基準で求めます。これを用いて

\[\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}\]

とすれば、実験と同じ基準の電極電位が得られます。

電気二重層容量

グランドポテンシャル\(\eqref{eq:omega-grand}\)を電極電位の関数として見ると、\(\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}\]

が成り立ちます。第1式は電極電荷 \(-e\Delta N\) についてのLippmannの式です。\(\Omega\) はPZCで極大をとる上に凸の放物線に近い形になり、その曲率から電気二重層容量 \(C_\mathrm{dl}\) が得られます。

適用例

Pt(111)/1 M HCl水溶液界面について、PZC、SHE電位、PZC近傍の電気二重層容量が実験値と整合することが示されています（Hagiwara et al. 2022）。再定式化したESM-RISM法では、PZCは \(+0.23\) V vs. SHEとなり、得られたSHE電位を絶対電位に換算すると4.06–4.37 Vで、IUPACの推奨値4.44 Vと整合します（Hagiwara and Otani 2026）。このほか、溶液の電荷分布と静電ポテンシャルの分布、吸着種のグランドポテンシャルの電位依存性から、電位に依存した吸着サイトや反応の平衡電位を議論できます。

参考文献

• 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