1D-RISM

1D-RISM（one-dimensional reference interaction site model）は、均一でバルクな溶液の構造を、分子を構成する原子サイト間の動径相関関数として求める積分方程式理論です。分子の配向について平均をとるため、相関関数はサイト間距離\(r\)だけの関数になります。ESM-RISM法の流れの中では、1D-RISMは溶媒どうしの相関を表す溶媒感受率\(\chi_{\alpha\gamma}\)を与える役割を担い、この\(\chi_{\alpha\gamma}\)が3D-RISMとESM-RISM法の入力になります。本章では原子単位（\(\hbar = e = m_\mathrm{e} = 1\)）を用い、\(\beta = 1/k_\mathrm{B}T\)とします。

分子性液体のOrnstein–Zernike方程式

単純液体の二体相関は、全相関関数\(h(r) = g(r) - 1\)と直接相関関数\(c(r)\)を結ぶOrnstein–Zernike（OZ）方程式で記述されます。

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

ここで\(\rho\)は数密度、\(g(r)\)は動径分布関数です。\(h\)は粒子1と2の間の相関を、直接の相関\(c\)と、第三の粒子を経由する間接的な相関の連鎖に分解したものです。Fourier空間では畳み込みが積になり、\(h(g) = c(g)/[1 - \rho c(g)]\)となります。

分子性液体では、位置\(\boldsymbol{r}_i\)と配向\(\Omega_i\)をまとめた座標\(i = (\boldsymbol{r}_i, \Omega_i)\)を用いて、分子OZ方程式が書けます。

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

\(\mathrm{d}(3) = \mathrm{d}\boldsymbol{r}_3\,\mathrm{d}\Omega_3\)で、非線形分子では\(\int \mathrm{d}\Omega = 8\pi^2\)です。この方程式は相関関数が6次元の変数に依存するため、そのまま数値的に解くのは困難です。

サイト–サイト表現：RISM方程式

ChandlerとAndersenは、分子を相互作用サイト（通常は原子）の集まりとみなし、分子の配向について平均をとったサイト–サイト相関関数\(h_{\alpha\gamma}(r)\)、\(c_{\alpha\gamma}(r)\)で液体の構造を記述しました。ここでギリシャ文字\(\alpha, \gamma\)は溶液を構成するすべての成分（溶媒分子とイオン）のサイトを表します。HirataとRosskyは、これを電荷をもつサイトからなる極性液体に拡張しました（extended RISM、XRISM）。行列形式で書いたRISM方程式は次のとおりです。

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

\(h\)、\(c\)、\(\omega\)はサイト添字をもつ行列、\(\rho\)は対角要素\(\rho_\alpha\)（サイト\(\alpha\)をもつ成分の数密度）からなる対角行列、\(*\)は畳み込み積分を表します。畳み込みはFourier変換によりFourier空間での行列積として評価します。

分子内相関関数

\(\omega_{\alpha\gamma}\)は分子内相関関数で、分子の幾何構造を表します。サイト\(\alpha\)と\(\gamma\)が同じ分子に属し、その間の距離が\(l_{\alpha\gamma}\)に固定されているとき、Fourier空間では

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

となります。これは実空間の\(\omega_{\alpha\gamma}(r) = \delta(r - l_{\alpha\gamma})/(4\pi l_{\alpha\gamma}^2)\)、すなわち半径\(l_{\alpha\gamma}\)の球殻に対応し、\(\alpha = \gamma\)（\(l_{\alpha\alpha} = 0\)）では\(\omega_{\alpha\alpha}(g) = 1\)です。\(\alpha\)と\(\gamma\)が異なる分子に属するときは\(\omega_{\alpha\gamma} = 0\)です。単原子イオンは\(\omega_{\alpha\alpha} = 1\)だけをもつ1サイト分子として扱います。

溶媒感受率

RISM方程式\(\eqref{eq:rism}\)の第1式の右辺は\(\omega * c * (\omega + \rho h)\)とまとめられます。括弧内の量が溶媒感受率です。

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

\(\chi_{\alpha\gamma}\)は、サイト\(\gamma\)に加えた摂動に対してサイト\(\alpha\)の密度がどう応答するかを、分子内（\(\omega\)）と分子間（\(\rho h\)）の両方の寄与を含めて表しています。収束した1D-RISMの解から得られるこの関数が、3D-RISMとESM-RISMで用いる量です。

相互作用ポテンシャルと閉包関係式

RISM方程式は\(h\)と\(c\)の2つの未知関数を含むため、もう1つの関係式（閉包関係式）が必要です。厳密な閉包関係式は

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

と書けます。\(g_{\alpha\gamma} = h_{\alpha\gamma} + 1\)は対分布関数、\(b_{\alpha\gamma}\)は高次の相関を含むブリッジ関数です。ブリッジ関数は一般に未知なので、近似を導入します。

サイト間の対ポテンシャル\(u_{\alpha\gamma}\)は、12-6型のLennard-Jones（LJ）ポテンシャルとクーロンポテンシャルの和で与えます。

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

\(q_\alpha\)はサイト電荷です。異種サイト間のLJパラメータはLorentz–Berthelot則\(\sigma_{\alpha\gamma} = (\sigma_\alpha + \sigma_\gamma)/2\)、\(\varepsilon_{\alpha\gamma} = \sqrt{\varepsilon_\alpha \varepsilon_\gamma}\)で決めます。

HNC閉包

ブリッジ関数を無視した近似がhypernetted-chain（HNC）閉包です。

\[\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閉包は静電相互作用の強い系にもよく用いられますが、引力の強いサイトの近くでは指数関数が大きくなりすぎ、反復計算が発散しやすくなります。

KH閉包

KovalenkoとHirataは、指数の引数が正になる領域だけHNC閉包を線形化した閉包（KH閉包）を提案しました。

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

分布関数が減少する領域（\(d \le 0\)）ではHNCと同じで、分布が密集する領域では線形に抑えられます。このため、金属表面や電荷をもつ溶質の近傍でも安定に解が得られます。ESM-RISM法の計算では、主にKH閉包を用います。

長距離クーロン相互作用の分離

直接相関関数は遠方で\(c_{\alpha\gamma}(r) \to -\beta u_{\alpha\gamma}(r)\)とふるまうため、クーロン相互作用に由来する\(1/r\)の裾をもちます。この裾をそのまま数値的にFourier変換すると精度と安定性が損なわれるので、\(c\)を短距離部分と長距離部分に分けます。

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

長距離部分は、幅\(\tau\)のGauss関数で遮蔽したクーロンポテンシャルで与えます。

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

\(u^\mathrm{L}\)は解析的に与えられるので反復の前に一度だけ計算し、反復では短距離部分\(c^\mathrm{S}\)だけを最適化します。\(u^\mathrm{S} = u - u^\mathrm{L}\)とおくと\(-\beta u + h - c = -\beta u^\mathrm{S} + h - c^\mathrm{S}\)が成り立つので、閉包関係式\(\eqref{eq:hnc}\)、\(\eqref{eq:kh}\)は短距離量だけで評価できます。RISM方程式\(\eqref{eq:rism}\)は、\(c^\mathrm{L}\)の解析的なFourier成分を加えてFourier空間で解きます。\(\tau\)としては1 Å程度が典型的に用いられます。

誘電的に整合したRISM（DRISM）

XRISMは、双極子をもつ分子からなる溶液の比誘電率\(\varepsilon\)を系統的に過小評価することが知られています。この欠点は、双極子–双極子相関のような高次の相関を担うブリッジ関数を無視したことに由来します。PerkynsとPettittが提案したDRISM（dielectrically consistent RISM）は、\(\varepsilon\)を経験的パラメータとして与え、溶液中の電荷間のクーロン相互作用が正しく遮蔽されるようにRISM方程式を修正します。

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

\(\xi_{\alpha\gamma}\)は誘電補正関数で、\(\alpha\)と\(\gamma\)がともに極性分子のサイトであるときだけ値をもち、イオンが関わる組では0です。

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

\(j_0\)、\(j_1\)は0次と1次の球Bessel関数、\(d^i_\alpha\)（\(i = x, y, z\)）は、双極子を\(z\)軸に向けて置いた極性分子の中で、電荷の絶対値の中心から測ったサイト\(\alpha\)の座標です。核関数\(\tilde{h}(g)\)は

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

です。\(\rho_\mathrm{tot}\)は溶液の全数密度、\(\rho_\mathrm{pol}\)は極性分子の数密度、\(\mu\)は分子1個あたりの双極子モーメント、\(y\)は双極子密度です。\(a\)は補正を長波長に限るためのパラメータで、水溶液では0.5 Å程度を用います。\(\varepsilon - 1 = 3y\)となる相関のない極限では補正が消えます。

\(\omega' + \rho h' = \omega + \rho h\)が成り立つので、DRISMの解からも式\(\eqref{eq:chi}\)をそのまま使って\(\chi_{\alpha\gamma}\)が得られます。つまり、誘電補正の効果は溶媒感受率を通して、3D-RISMとESM-RISMにそのまま引き継がれます。HagiwaraらはDRISMをESM-RISMと組み合わせ、NaCl水溶液の平均活量係数、Pt(111)/HCl水溶液界面のゼロ電荷電位や電気二重層容量がXRISMより改善することを示しました。

過剰化学ポテンシャル

1D-RISMの解から、溶液中の成分\(s\)の過剰化学ポテンシャル（溶媒和自由エネルギー）を閉形式で計算できます。KH閉包に対しては次式になります。

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

\(\Theta\)はHeavisideの階段関数です。HNC閉包では\(\Theta\)を外した式（Singer–Chandlerの式）に、Gauss揺らぎ（GF）近似では\(\frac{1}{2}h^2\)の項を落とした式になります。この表式は閉包関係式と整合したもので、熱力学積分を行わずに一度の計算で自由エネルギーが求まる点がRISM理論の利点です。

有限濃度の電解質溶液では、平均活量係数\(\gamma_{m,\pm}\)を過剰化学ポテンシャルから求められます。

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

\(i\)はイオン種、\(s_i\)は化学量論係数で\(\nu = \sum_i s_i\)、\(\rho_\mathrm{w}\)と\(\rho_{\mathrm{w},\infty}\)は有限濃度と無限希釈の溶液における水の数密度です。XRISMは\(\gamma_{m,\pm}\)を大きく過小評価しますが、DRISMでは実験値をよく再現します。

数値解法

1D-RISMの相関関数は動径関数なので、Fourier変換は球Bessel変換（サイン変換）になります。

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

計算は次の反復で進めます。

1. \(c^\mathrm{S}_{\alpha\gamma}(r)\)の初期値を与え、Fourier変換して\(c^\mathrm{L}\)を加える。
2. Fourier空間でRISM方程式\(\eqref{eq:rism}\)（またはDRISMの式\(\eqref{eq:drism}\)）を各\(g\)ごとの行列演算で解き、\(h_{\alpha\gamma}(g)\)を得る。
3. 逆Fourier変換した\(h_{\alpha\gamma}(r)\)を閉包関係式に代入し、新しい\(c^\mathrm{S}_{\alpha\gamma}(r)\)を得る。
4. 残差が収束判定値を下回るまで繰り返す。

単純な混合による反復は収束が遅いため、modified direct inversion in the iterative subspace（MDIIS）法で加速します。Hagiwaraらの計算では、実空間カットオフ999.8 Bohr、逆空間カットオフ18.763 Bohr\(^{-1}\)、5973点の格子を用い、MDIISの収束判定値を\(10^{-8}\)としています。

3D-RISMとESM-RISMへの受け渡し

1D-RISMの主な出力は、収束した\(h_{\alpha\gamma}(g)\)から式\(\eqref{eq:chi}\)で組み立てた溶媒感受率\(\chi_{\alpha\gamma}(g)\)です。これはバルク溶液の性質だけで決まり、溶質には依存しません。

• 3D-RISMでは、\(\chi_{\alpha\gamma}(|\boldsymbol{g}|)\)を3次元の逆格子ベクトル上に補間して、溶質まわりの3次元分布を求めます。
• ESM-RISM法では、\(\chi_{\alpha\gamma}(g)\)を面直方向\(g_z\)について1次元逆Fourier変換し、Laue表示\(\chi_{\alpha\gamma}(g_\parallel, |z - z'|)\)にして用います。\(\omega_{\alpha\gamma}(g)\)は\(g\)の大きいところまで減衰しないので、この変換には工夫が必要です（詳細はESM-RISM法の章を参照してください）。

溶媒のモデル（サイト構造、LJパラメータ、電荷）、温度、濃度、誘電補正の有無はすべてこの段階で決まるため、1D-RISMの精度が3D-RISMとESM-RISMの結果を左右します。

参考文献

• 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