コンテンツにスキップ

3D-RISM

3D-RISM（three-dimensional reference interaction site model）は、溶質のまわりの溶媒サイトの分布を3次元の関数として求める積分方程式理論です。1D-RISMが溶質を配向平均した動径関数で扱うのに対し、3D-RISMでは溶質を空間に固定し、その形状と電荷分布を3次元のまま取り込みます。溶質を密度汎関数理論（DFT）で扱う電子・イオン系とすれば、DFTと溶液理論を自己無撞着に結合したハイブリッド溶媒和モデルになります。本章では原子単位（\(\hbar = e = m_\mathrm{e} = 1\)）を用い、\(\beta = 1/k_\mathrm{B}T\)とします。溶媒サイトを\(\alpha, \gamma\)、溶質原子を\(I\)で表します。

1D-RISMから3D-RISMへ

溶質（添字\(a, b\)）と溶媒（添字\(\alpha, \gamma, \lambda\)）からなる混合系を考え、溶質を無限希釈とします。サイト–サイトRISM方程式\(h = \omega * c * \chi\)の溶質–溶媒成分は

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

と書けます。ここで\(\mu, \lambda\)は溶質と溶媒のすべてのサイトを走ります。3D-RISMでは、溶質の構造を分子内相関関数による配向平均で表すのではなく、溶質を空間に固定してその構造を3次元の相関関数に直接もたせます。そこで溶質サイトについて

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

とおき、溶質サイトどうしの直接相関を0とします。すると\(\mu\)の和は\(\mu = a\)だけになり、\(\lambda\)の和は溶媒サイトに限られます。溶質全体をひとつの3次元の対象とみなして添字\(a\)を省き、溶媒どうしの感受率\(\chi_{\alpha\gamma}\)を溶質に乱されないバルクの等方的な関数（1D-RISMの解）で近似すると、3D-RISM方程式が得られます。

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

\(h_\gamma(\boldsymbol{r}) = g_\gamma(\boldsymbol{r}) - 1\)と\(c_\gamma(\boldsymbol{r})\)は、溶質に対する溶媒サイト\(\gamma\)の3次元の全相関関数と直接相関関数です。\(g_\gamma(\boldsymbol{r})\)は位置\(\boldsymbol{r}\)におけるサイト\(\gamma\)の密度をバルク密度\(\rho_\gamma\)で規格化したものです。\(\chi_{\alpha\gamma}(g) = \omega_{\alpha\gamma}(g) + \rho_\alpha h_{\alpha\gamma}(g)\)は1D-RISMで求めた溶媒感受率です。畳み込みはFourier空間で積になります。

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

溶質の情報は\(c_\alpha\)と次に述べる相互作用ポテンシャル\(u_\gamma\)だけに入り、溶媒の情報は\(\chi_{\alpha\gamma}\)だけに入ります。この分離のおかげで、溶媒の計算（1D-RISM）を一度行えば、さまざまな溶質に使い回せます。

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

溶媒サイト\(\gamma\)が位置\(\boldsymbol{r}\)で感じるポテンシャル\(u_\gamma(\boldsymbol{r})\)は、溶質原子とのLennard-Jones（LJ）相互作用と、溶質の電荷分布がつくる静電ポテンシャルの和です。

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

\(\boldsymbol{R}_I\)は溶質原子\(I\)の位置、\(q_\gamma\)は溶媒サイトの電荷です。LJパラメータはLorentz–Berthelot則で組み合わせます。\(v_\mathrm{DFT}\)はDFTの静電ポテンシャルで、Hartreeポテンシャルと擬ポテンシャルの局所部分の和です。

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

\(n(\boldsymbol{r})\)は電子密度です。\(v_\mathrm{DFT}\)は電子（電荷\(-1\)）に対するポテンシャルエネルギーなので、電荷\(q_\gamma\)のサイトが感じる静電エネルギーは\(-q_\gamma v_\mathrm{DFT}\)になります。電子の分布が変わると\(u_\gamma\)が変わり、溶媒の分布が変わります。これがDFTから溶媒への結合です。

閉包関係式

3D-RISM方程式も閉包関係式と連立して解きます。HNC閉包は

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

です。金属電極のように溶媒を強く引きつける溶質ではHNCの指数関数が大きくなりすぎるため、KovalenkoとHirataが金属–水界面の計算のために導入したKH閉包を用います。

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

溶質内部ではLJの斥力により\(g_\gamma \to 0\)となり、溶媒は溶質の外側にだけ分布します。

長距離漸近形と電荷の扱い

溶質から遠いところでは\(c_\gamma(\boldsymbol{r}) \to -\beta u_\gamma(\boldsymbol{r})\)となるので、溶質が電荷や双極子をもつと\(c_\gamma\)はクーロン型の長い裾をもちます。1D-RISMと同様に、\(c_\gamma\)を短距離部分と長距離部分に分けます。

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

長距離部分は、DFTの静電ポテンシャルを幅\(\tau\)のGauss関数でぼかしたものとしてFourier空間で与えます。

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

\(\tau\)としては1 Å程度が典型的です。\(u^\mathrm{L}_\gamma\)は電子密度が与えられれば3D-RISMの反復の前に決まるので、反復では短距離部分\(c^\mathrm{S}_\gamma\)だけを最適化します。閉包関係式の指数は\(-\beta u^\mathrm{S}_\gamma + h_\gamma - c^\mathrm{S}_\gamma\)（\(u^\mathrm{S}_\gamma = u_\gamma - u^\mathrm{L}_\gamma\)）と短距離量だけで書け、3D-RISM方程式\(\eqref{eq:3drism-g}\)には\(c^\mathrm{S}_\alpha + c^\mathrm{L}_\alpha\)を代入します。

電荷をもつ溶質では、全相関関数\(h_\gamma\)も遠方でDebye型に遮蔽された長い裾をもちます。孤立した溶質を扱う3D-RISMでは、\(h\)と\(c\)の両方の漸近形を解析的に分離して扱う方法が開発されています。周期的スーパーセルでは\(\boldsymbol{g} = 0\)のクーロン成分が定まらないため、通常の平面波DFTと同様に、静電ポテンシャルの平均値を0にとり、系に正味の電荷があれば一様な背景電荷で中和する必要があります。

溶媒和自由エネルギー

3D-RISMの解から、溶媒和自由エネルギー（過剰化学ポテンシャル）\(\Delta\mu_\mathrm{solv}\)を閉形式で計算できます。KH閉包に対応する表式は次のとおりです。

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

HNC閉包では\(\Theta\)を外した式になります。Gauss揺らぎ（GF）近似では\(\frac{1}{2}h^2\)の項を落とします。

\[\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らの計算では、KH閉包で分布を求め、自由エネルギーはGF表式で評価しています。KH閉包とHNC閉包の表式は、解において\(h_\gamma\)と\(c_\gamma\)について停留的です。この性質は、次節で電子密度に関する汎関数微分を求めるときに使います。

DFTとの結合

Helmholtz自由エネルギー

DFT-RISM法では、DFTの全エネルギー\(E_\mathrm{DFT}\)を、溶媒和自由エネルギーを加えたHelmholtz自由エネルギー\(A\)に置き換えます。

\[\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]\)は非相互作用運動エネルギー、外部ポテンシャルのエネルギー、Hartreeエネルギー、交換相関エネルギーの和です。溶質の原子配置と電子数を固定して\(A\)を電子密度について変分すると、溶媒中の溶質の電子状態が得られます。

Kohn–Shamポテンシャルへのフィードバック

Kohn–Sham有効ポテンシャルには、溶媒からの寄与\(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}\]

自由エネルギー\(\eqref{eq:mu-kh3d}\)が\(h_\gamma\)、\(c_\gamma\)について停留的であることから、\(\delta\Delta\mu_\mathrm{solv}/\delta u_\gamma(\boldsymbol{r}) = \rho_\gamma g_\gamma(\boldsymbol{r})\)が成り立ちます。式\(\eqref{eq:u3d}\)、\(\eqref{eq:vdft}\)から\(\delta u_\gamma(\boldsymbol{r}')/\delta n(\boldsymbol{r}) = -q_\gamma/|\boldsymbol{r} - \boldsymbol{r}'|\)なので、

\[\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}\)は溶媒の電荷密度で、\(v_\mathrm{solv}\)はそれがつくる静電ポテンシャルを電子のポテンシャルエネルギーとして表したものです。LJ相互作用は電子密度に依存しないので、\(v_\mathrm{solv}\)には現れません。計算は次の自己無撞着ループで進めます。

  1. 電子密度\(n\)から\(v_\mathrm{DFT}\)を求め、式\(\eqref{eq:u3d}\)で\(u_\gamma\)をつくる。
  2. 3D-RISM方程式と閉包関係式を解き、\(g_\gamma\)を得る。
  3. 式\(\eqref{eq:vsolv}\)で\(\rho_\mathrm{solv}\)と\(v_\mathrm{solv}\)を求め、Kohn–Sham方程式を解いて新しい\(n\)を得る。
  4. \(n\)と\(g_\gamma\)がともに収束するまで繰り返す。

Janakの定理

Kohn–Sham軌道の占有数を\(f_i\)、固有値を\(\epsilon_i\)とすると、Helmholtz自由エネルギー\(\eqref{eq:helmholtz}\)についてもJanakの定理が成り立ちます。

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

証明は通常のDFTと同じ手順で、\(\partial \Delta\mu_\mathrm{solv}/\partial f_i = \int \mathrm{d}\boldsymbol{r}\, |\psi_i(\boldsymbol{r})|^2 v_\mathrm{solv}(\boldsymbol{r})\)が\(v_\mathrm{eff}\)の一部として吸収されることによります。したがって電子の化学ポテンシャルは\(\mu_\mathrm{e} = \partial A/\partial N\)で与えられ、溶媒中の電極電位を電子の化学ポテンシャルから評価する根拠になります。

周期的スーパーセルでの解法

平面波基底のDFTと組み合わせる場合、3D-RISM方程式はDFTと同じ周期的スーパーセルで解きます。

  • 相関関数は実空間の格子と逆格子ベクトル\(\boldsymbol{g}\)の間をFFTで行き来します。Hagiwaraらの計算では、相関関数のカットオフエネルギーを160 Ryとしています。
  • 1D-RISMで動径格子上に求めた\(\chi_{\alpha\gamma}(g)\)を、各逆格子ベクトルの長さ\(|\boldsymbol{g}|\)に補間して式\(\eqref{eq:3drism-g}\)に用います。
  • 閉包関係式は実空間の各格子点で、3D-RISM方程式は逆空間の各\(\boldsymbol{g}\)で評価し、MDIIS法で\(c^\mathrm{S}_\gamma\)の反復を加速します。

この定式化はカノニカル条件に基づきます。セル内の溶媒サイト\(\gamma\)の数

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

は計算の間一定に保たれます。この性質は、層状物質の層間に入る水分子の数を評価するといった用途に役立ちます。

界面への適用と限界

3D-RISMは溶液中の分子やイオンには有効ですが、周期的スーパーセルのままでは電極/電解液界面を扱うのが困難です。

  • スラブとその間の溶液が面直方向にも周期的に繰り返されるため、溶液領域は有限の厚さしかもたず、バルク溶液（溜め）につながっていません。電位の基準となるバルク溶液の内部電位を定義できません。
  • 溶媒粒子の数が一定のカノニカル条件なので、電極の帯電に応じて溶液側から対イオンが出入りする電気二重層の形成を記述できません。
  • 帯電した電極を扱うには一様な背景電荷が必要になり、界面の静電ポテンシャルとエネルギーに人為的な影響が入ります。

ESM-RISM法はこれらの問題を解決します。3D-RISM方程式を面内方向だけFourier変換したLaue表示に書き直し、

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

面直方向には溶液が無限遠のバルクまで続く境界条件を課します。静電ポテンシャルはESM法で開放境界条件のもとで解きます。これにより電子と溶媒粒子の両方をグランドカノニカル条件で扱え、グランドポテンシャル

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

を用いて、背景電荷なしに電極電位を制御した計算ができます。ここで\(N_0\)は中性系の電子数です。

参考文献

  • 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