ESM法¶

有効遮蔽媒質（effective screening medium: ESM）法は、帯電した表面・界面や電圧を印加した表面を、平面波基底の第一原理計算で扱うための静電ポテンシャルの取り扱いです。スラブに垂直な方向（\(z\)方向）の周期境界条件を外し、スラブの外側に置いた連続媒質（真空・理想金属など）が与える境界条件の下でPoisson方程式をグリーン関数で解きます。Kohn-Sham方程式は通常の周期境界条件の下で解くため、既存の平面波・擬ポテンシャルコードへの変更は静電相互作用の部分に限られます。電極電位を一定に保つ計算はConstant-μ法で扱います。

帯電表面と周期境界条件¶

1990年代中頃から、燃料電池を主な対象として、電圧を印加された電極／電解質界面（電気化学界面）を第一原理計算で扱う必要性が高まりました。界面には電気二重層が形成され、電極表面の過剰電荷とそれを打ち消す対電荷の間に電場が生じます。

平面波基底の計算では、計算セルを3次元的に繰り返す周期境界条件を課します。この条件の下で電気二重層を表現しようとすると、セルの端で静電ポテンシャルが不連続になります。セル全体が帯電していれば静電エネルギーは発散し、これを避けるために一様背景電荷を導入すると、電子のエネルギー準位の絶対値が物理的な意味を失います。鋸型ポテンシャルや一様背景電荷を用いる従来の方法は、電気化学界面のモデルとしては不十分でした。

ESM法では次のように扱います。

表面平行方向（\(xy\)方向）の周期境界条件は維持し、\(z\)方向の周期境界条件を外す。
スラブの外側に誘電率\(\varepsilon(\boldsymbol{r})\)で特徴づけられる遮蔽媒質を置き、モデルに応じた境界条件を課す（孤立スラブ、電圧を印加した表面、キャパシターなど）。
遮蔽媒質がスラブの過剰電荷を打ち消す電荷を担うため、媒質まで含めた系全体は常に電気的に中性になる。

全エネルギー汎関数と一般化Poisson方程式¶

電子の電荷密度を\(\rho_\mathrm{e}\)、イオンの電荷密度を\(\rho_\mathrm{I}\)、全電荷密度を\(\rho_\mathrm{tot}=\rho_\mathrm{e}+\rho_\mathrm{I}\)とします。平面波・擬ポテンシャル法の慣例に合わせて、電子の電荷を正にとります。静電ポテンシャル\(V(\boldsymbol{r})\)も独立変数とみなし、全エネルギーを次の汎関数で定義します。

\[\begin{equation} E[\rho_\mathrm{e},V]=T_\mathrm{s}[\rho_\mathrm{e}]+E_\mathrm{xc}[\rho_\mathrm{e}]-\int d\boldsymbol{r}\,\frac{\varepsilon(\boldsymbol{r})}{8\pi}\left|\nabla V(\boldsymbol{r})\right|^2+\int d\boldsymbol{r}\,\rho_\mathrm{tot}(\boldsymbol{r})V(\boldsymbol{r}) \label{eq:functional} \end{equation}\]

ここで\(T_\mathrm{s}\)は運動エネルギー、\(E_\mathrm{xc}\)は交換相関エネルギー、\(\varepsilon(\boldsymbol{r})\)は位置に依存する比誘電率です。\(V\)について変分をとると

\[\begin{equation} \delta E=\int d\boldsymbol{r}\left[-\frac{\varepsilon(\boldsymbol{r})}{4\pi}\nabla V\cdot\nabla\delta V+\rho_\mathrm{tot}\,\delta V\right] =-\oint_S\frac{\varepsilon(\boldsymbol{r})}{4\pi}\,\delta V\,\nabla V\cdot d\boldsymbol{S}+\int d\boldsymbol{r}\left[\nabla\cdot\left(\frac{\varepsilon(\boldsymbol{r})}{4\pi}\nabla V\right)+\rho_\mathrm{tot}\right]\delta V \label{eq:variation} \end{equation}\]

となります。境界条件により表面積分は消えるので、停留条件\(\delta E/\delta V=0\)から一般化Poisson方程式と、対応するグリーン関数の方程式が得られます。

\[\begin{align} \nabla\cdot\left[\varepsilon(\boldsymbol{r})\nabla\right]V(\boldsymbol{r})&=-4\pi\rho_\mathrm{tot}(\boldsymbol{r}),\label{eq:poisson}\\ \nabla\cdot\left[\varepsilon(\boldsymbol{r})\nabla\right]G(\boldsymbol{r},\boldsymbol{r}')&=-4\pi\delta(\boldsymbol{r}-\boldsymbol{r}').\label{eq:green} \end{align}\]

静電ポテンシャルはグリーン関数を用いて

\[\begin{equation} V(\boldsymbol{r})=\int d\boldsymbol{r}'\,G(\boldsymbol{r},\boldsymbol{r}')\rho_\mathrm{tot}(\boldsymbol{r}') \label{eq:v-green} \end{equation}\]

と表せます。Poisson方程式の解を式\(\eqref{eq:functional}\)に代入し、恒等式\(\nabla\cdot(\varepsilon V\nabla V)=V\nabla\cdot(\varepsilon\nabla V)+\varepsilon|\nabla V|^2\)を用いて部分積分すると、静電項は\(\frac{1}{2}\int d\boldsymbol{r}\,\rho_\mathrm{tot}V\)となり、全エネルギーは通常の形に戻ります。

\[\begin{equation} E[\rho_\mathrm{e}]=T_\mathrm{s}[\rho_\mathrm{e}]+E_\mathrm{xc}[\rho_\mathrm{e}]+\frac{1}{2}\iint d\boldsymbol{r}\,d\boldsymbol{r}'\,\rho_\mathrm{tot}(\boldsymbol{r})G(\boldsymbol{r},\boldsymbol{r}')\rho_\mathrm{tot}(\boldsymbol{r}') \label{eq:energy-green} \end{equation}\]

静電項は電子間（Hartree）、電子–イオン間、イオン間の3項に分かれます。

\[\begin{equation} E_\mathrm{H}=\frac{1}{2}\iint d\boldsymbol{r}\,d\boldsymbol{r}'\,\rho_\mathrm{e}(\boldsymbol{r})G(\boldsymbol{r},\boldsymbol{r}')\rho_\mathrm{e}(\boldsymbol{r}'),\quad E_\mathrm{e\text{-}I}=\iint d\boldsymbol{r}\,d\boldsymbol{r}'\,\rho_\mathrm{e}(\boldsymbol{r})G(\boldsymbol{r},\boldsymbol{r}')\rho_\mathrm{I}(\boldsymbol{r}'),\quad E_\mathrm{I\text{-}I}=\frac{1}{2}\iint d\boldsymbol{r}\,d\boldsymbol{r}'\,\rho_\mathrm{I}(\boldsymbol{r})G(\boldsymbol{r},\boldsymbol{r}')\rho_\mathrm{I}(\boldsymbol{r}') \label{eq:energy-terms} \end{equation}\]

つまりESM法は、DFTの全エネルギーの形を保ったまま、静電相互作用を\(1/|\boldsymbol{r}-\boldsymbol{r}'|\)からグリーン関数\(G(\boldsymbol{r},\boldsymbol{r}')\)に置き換える方法です。

式\(\eqref{eq:functional}\)を軌道について変分すると、Kohn-Sham方程式が得られます。

\[\begin{equation} \left[-\frac{1}{2}\nabla^2+V(\boldsymbol{r})+\hat{V}_\mathrm{NL}+V_\mathrm{xc}(\boldsymbol{r})\right]\psi_i(\boldsymbol{r})=\epsilon_i\psi_i(\boldsymbol{r}) \label{eq:ks} \end{equation}\]

2つのEuler-Lagrange方程式は別々に扱います。

Poisson方程式：モデルに応じた境界条件の下で、あらかじめ解析的に求めたグリーン関数を用いて解きます。グリーン関数はプログラムに直接書き込むので、計算の実行中に微分方程式を解く必要はなく、計算負荷の増加は最小限に抑えられます。
Kohn-Sham方程式：通常の周期境界条件のソルバーをそのまま使います。得られた静電ポテンシャルは\(z\)方向に周期的ではないため、電子密度が無視できるセルの端の近くで滑らかにつなぎます。この扱いは、スラブの電子密度がセルの端に達するまでに十分減衰していれば正当です。

擬ポテンシャルは真空中の孤立原子に対して作られているため、原子核は\(\varepsilon=1\)の領域に置きます。

Laue表示¶

\(\varepsilon\)が\(z\)のみに依存する場合を考えます。\(xy\)方向には周期境界条件が課されているので、Poisson方程式は面内方向にFourier変換したLaue表示で扱えます。

\[\begin{equation} V(\boldsymbol{g}_\parallel,z)=\int d\boldsymbol{r}_\parallel\,V(\boldsymbol{r}_\parallel,z)\,e^{-i\boldsymbol{g}_\parallel\cdot\boldsymbol{r}_\parallel},\qquad V(\boldsymbol{r}_\parallel,z)=\frac{1}{(2\pi)^2}\int d\boldsymbol{g}_\parallel\,V(\boldsymbol{g}_\parallel,z)\,e^{i\boldsymbol{g}_\parallel\cdot\boldsymbol{r}_\parallel} \label{eq:laue-ft} \end{equation}\]

電荷密度とグリーン関数も同様に変換すると、Laue表示のPoisson方程式が得られます（\(g_\parallel=|\boldsymbol{g}_\parallel|\)）。

\[\begin{align} \left\{\partial_z\left[\varepsilon(z)\partial_z\right]-\varepsilon(z)g_\parallel^2\right\}V(\boldsymbol{g}_\parallel,z)&=-4\pi\rho_\mathrm{tot}(\boldsymbol{g}_\parallel,z),\label{eq:poisson-laue}\\ \left\{\partial_z\left[\varepsilon(z)\partial_z\right]-\varepsilon(z)g_\parallel^2\right\}G(g_\parallel,z,z')&=-4\pi\delta(z-z').\label{eq:green-laue} \end{align}\]

静電ポテンシャルは

\[\begin{equation} V(\boldsymbol{g}_\parallel,z)=\int dz'\,G(g_\parallel,z,z')\rho_\mathrm{tot}(\boldsymbol{g}_\parallel,z') \label{eq:v-laue} \end{equation}\]

で与えられます。問題は、各\(g_\parallel\)について\(z\)に関する1次元の2階常微分方程式を解くことに帰着しました。これは適切な境界条件の下で必ず解けます。

以下では、スラブの電子が\(-z_0\le z\le z_0\)（長さ\(L=2z_0\)の計算セル）に閉じ込められ、遮蔽媒質の境界が\(z=\pm z_1\)（\(z_1\ge z_0\)）にあるとします。

グリーン関数の構成¶

非同次Sturm-Liouville方程式

\[\begin{equation} \hat{L}y(z)+f(z)=0,\qquad \hat{L}=\frac{d}{dz}\,p(z)\frac{d}{dz}+q(z)\qquad (a\le z\le b) \label{eq:sturm} \end{equation}\]

のグリーン関数は、次の条件で決まります。

\(z\ne z'\)では同次方程式\(\hat{L}G=0\)を満たす。
\(z=a\)および\(z=b\)で、\(y\)に課す境界条件を満たす。
\(z=z'\)で連続である。
導関数が\(z=z'\)で\(\left.\partial_z G\right|_{z'+0}-\left.\partial_z G\right|_{z'-0}=-1/p(z')\)だけ跳ぶ。

左端の境界条件を満たす同次解を\(\psi_\mathrm{L}\)、右端の境界条件を満たす同次解を\(\psi_\mathrm{R}\)とします。Wronskianについて\(p(z)\left[\psi_\mathrm{L}\psi_\mathrm{R}'-\psi_\mathrm{L}'\psi_\mathrm{R}\right]=A\)は定数となり、条件3と4から

\[\begin{equation} G(z,z')=-\frac{1}{A}\,\psi_\mathrm{L}(z_<)\,\psi_\mathrm{R}(z_>),\qquad z_<=\min(z,z'),\quad z_>=\max(z,z') \label{eq:green-sl} \end{equation}\]

が得られます。この方法を用いれば、任意の境界条件に対してグリーン関数を構成できます。

\(\varepsilon=1\)の領域では、式\(\eqref{eq:poisson-laue}\)は

\[\begin{equation} \frac{1}{4\pi}\left(\partial_z^2-g_\parallel^2\right)V(\boldsymbol{g}_\parallel,z)+\rho_\mathrm{tot}(\boldsymbol{g}_\parallel,z)=0 \label{eq:poisson-vac} \end{equation}\]

となり、\(p=1/4\pi\)、\(q=-g_\parallel^2/4\pi\)、\(f=\rho_\mathrm{tot}\)とおけば式\(\eqref{eq:sturm}\)と同じ形です。同次解は\(e^{\pm g_\parallel z}\)の線形結合なので、各境界条件の違いは\(\psi_\mathrm{L}\)と\(\psi_\mathrm{R}\)の選び方だけに現れます。比誘電率が無限大の領域は理想金属（電極）を表し、その内部では電位が一定です。この電位をエネルギーの基準（ゼロ）にとります。

境界条件とグリーン関数¶

ESM法の代表的な境界条件は次の3つです。

境界条件	構成	主な用途
bc1	真空／スラブ／真空	中性表面、極性表面、仕事関数
bc2	金属／スラブ／金属	キャパシター中のスラブ、電極間の電圧印加
bc3	真空／スラブ／金属	STM探針・ゲート電極、対極をもつ電気化学界面

bc1：真空／スラブ／真空¶

境界条件は

\[\begin{equation} \left.\partial_z V(\boldsymbol{g}_\parallel,z)\right|_{z=\pm\infty}=0,\qquad \varepsilon(z)=1 \label{eq:bc1} \end{equation}\]

です。同次解\(\psi_\mathrm{L}=a\,e^{g_\parallel z}\)、\(\psi_\mathrm{R}=b\,e^{-g_\parallel z}\)から\(A=-2abg_\parallel/4\pi\)となり、

\[\begin{equation} G^\mathrm{bc1}(g_\parallel,z,z')=\frac{4\pi}{2g_\parallel}\,e^{-g_\parallel|z-z'|}\qquad(g_\parallel\ne0) \label{eq:g-bc1} \end{equation}\]

が得られます。これは裸のCoulomb相互作用\(1/|\boldsymbol{r}-\boldsymbol{r}'|\)のLaue表示そのものです。電荷を遮蔽するものがないため、\(g_\parallel=0\)では発散します。したがってbc1は電気的に中性な系にしか使えません。帯電していると\(\partial_zV(0,z)\)が\(z\to\pm\infty\)の少なくとも一方でゼロにならず、ポテンシャルが\(z\)に比例して増減し、静電エネルギーが発散します。

bc1では静電ポテンシャルが常に真空準位（無限遠）を基準とするため、仕事関数がそのまま求まります。スラブが双極子モーメントをもちポテンシャルに長い裾が生じる場合も、繰り返しスラブ計算のような双極子補正は不要です。

bc2：金属／スラブ／金属¶

両側の\(|z|\ge z_1\)に理想金属を置きます。電圧を印加しない場合の境界条件は

\[\begin{equation} V(\boldsymbol{g}_\parallel,\pm z_1)=0,\qquad \varepsilon(z)=\begin{cases}1 & (|z|<z_1)\\ \infty & (|z|\ge z_1)\end{cases} \label{eq:bc2} \end{equation}\]

です。同次解\(\psi_\mathrm{L}=a\sinh[g_\parallel(z+z_1)]\)、\(\psi_\mathrm{R}=b\sinh[g_\parallel(z-z_1)]\)から\(A=abg_\parallel\sinh(2g_\parallel z_1)/4\pi\)となり、

\[\begin{equation} G^\mathrm{bc2}(g_\parallel,z,z')=\frac{4\pi}{2g_\parallel}\,e^{-g_\parallel|z-z'|}+\frac{4\pi}{2g_\parallel}\,\frac{e^{-2g_\parallel z_1}\cosh[g_\parallel(z-z')]-\cosh[g_\parallel(z+z')]}{\sinh(2g_\parallel z_1)} \label{eq:g-bc2} \end{equation}\]

が得られます（\(|z|,|z'|<z_1\)）。第1項は裸のCoulomb項、第2項は両電極に誘起される鏡像電荷の寄与です。第2項を展開すると

\[\begin{equation} \frac{4\pi}{2g_\parallel}\sum_{m=0}^{\infty}\left\{e^{-g_\parallel[z-z'+(4m+4)z_1]}+e^{g_\parallel[z-z'-(4m+4)z_1]}-e^{-g_\parallel[z+z'+(4m+2)z_1]}-e^{g_\parallel[z+z'-(4m+2)z_1]}\right\} \label{eq:g-bc2-series} \end{equation}\]

となり、\(z=z'\pm(4m+4)z_1\)にある同符号の鏡像電荷と、\(z=-z'\pm(4m+2)z_1\)にある異符号の鏡像電荷の無限列を表すことがわかります。

bc3：真空／スラブ／金属¶

片側の\(z\ge z_1\)にだけ理想金属を置きます。境界条件は

\[\begin{equation} V(\boldsymbol{g}_\parallel,z_1)=0,\qquad \left.\partial_z V(\boldsymbol{g}_\parallel,z)\right|_{z=-\infty}=0,\qquad \varepsilon(z)=\begin{cases}1 & (z<z_1)\\ \infty & (z\ge z_1)\end{cases} \label{eq:bc3} \end{equation}\]

です。同次解\(\psi_\mathrm{L}=a\,e^{g_\parallel z}\)、\(\psi_\mathrm{R}=b\sinh[g_\parallel(z-z_1)]\)から\(A=abg_\parallel e^{g_\parallel z_1}/4\pi\)となり、

\[\begin{equation} G^\mathrm{bc3}(g_\parallel,z,z')=\frac{4\pi}{2g_\parallel}\,e^{-g_\parallel|z-z'|}-\frac{4\pi}{2g_\parallel}\,e^{-g_\parallel(2z_1-z-z')} \label{eq:g-bc3} \end{equation}\]

が得られます（\(z,z'<z_1\)）。bc2と異なり、鏡像電荷は\(z=2z_1-z'\)の1つだけです。この境界条件は、電極（スラブ）と対極（ESM金属）からなる電気化学セルのモデルであり、Constant-μ法と組み合わせて用います。

面内波数ゼロの成分¶

\(g_\parallel=0\)では式\(\eqref{eq:green-laue}\)は1次元Poisson方程式\(\partial_z^2G=-4\pi\delta(z-z')\)となり、解は\(-2\pi|z-z'|\)に1次関数を加えたものです。各境界条件を課すと

\[\begin{align} G^\mathrm{bc1}(0,z,z')&=-2\pi|z-z'|+\mathrm{const.},\label{eq:g0-bc1}\\ G^\mathrm{bc2}(0,z,z')&=2\pi\left(z_1-\frac{zz'}{z_1}-|z-z'|\right),\label{eq:g0-bc2}\\ G^\mathrm{bc3}(0,z,z')&=2\pi\left(2z_1-z-z'-|z-z'|\right)\label{eq:g0-bc3} \end{align}\]

となります。bc2とbc3の式は、それぞれ式\(\eqref{eq:g-bc2}\)と式\(\eqref{eq:g-bc3}\)の\(g_\parallel\to0\)の極限とも一致します。bc2とbc3では金属電極が対電荷を担うため、スラブが帯電していても\(g_\parallel=0\)成分は有限です。bc1の定数項は任意ですが、系が中性であれば全エネルギーに寄与しません。

実空間表示と鏡像電荷¶

面内の距離を\(R=|\boldsymbol{r}_\parallel-\boldsymbol{r}'_\parallel|\)とし、\(\frac{4\pi}{2g_\parallel}e^{-g_\parallel|\zeta|}\)が\(1/\sqrt{R^2+\zeta^2}\)の2次元Fourier変換であることを使うと、実空間のグリーン関数が得られます。

\[\begin{align} G^\mathrm{bc1}&=\frac{1}{\sqrt{R^2+(z-z')^2}},\label{eq:g-real-bc1}\\ G^\mathrm{bc2}&=\frac{1}{\sqrt{R^2+(z-z')^2}}+\sum_{m=0}^{\infty}\left[\frac{1}{\sqrt{R^2+(z-z'+(4m+4)z_1)^2}}+\frac{1}{\sqrt{R^2+(z-z'-(4m+4)z_1)^2}}\right.\nonumber\\ &\qquad\left.-\frac{1}{\sqrt{R^2+(z+z'+(4m+2)z_1)^2}}-\frac{1}{\sqrt{R^2+(z+z'-(4m+2)z_1)^2}}\right],\label{eq:g-real-bc2}\\ G^\mathrm{bc3}&=\frac{1}{\sqrt{R^2+(z-z')^2}}-\frac{1}{\sqrt{R^2+(2z_1-z-z')^2}}.\label{eq:g-real-bc3} \end{align}\]

いずれも「裸のCoulomb項＋鏡像電荷項」の形をしており、ESM法の静電相互作用が古典電磁気学の鏡像法と同じ内容をもつことがわかります。電気力線はESM金属の内部には入りません。

全エネルギーの各項¶

平面波・擬ポテンシャル法では、グリーン関数が現れる項、すなわちイオン間相互作用（Ewald和）、局所擬ポテンシャル、Hartreeポテンシャルを修正します。非局所擬ポテンシャル、交換相関、運動エネルギーは変更不要です。OpenMXやSIESTAのように全電荷密度から静電ポテンシャルを一度に求める構成のコードでは、その部分の修正だけで済みます。

Hartreeポテンシャル¶

電子密度を計算セル内で\(\rho_\mathrm{e}(\boldsymbol{g}_\parallel,z)=\sum_{g_z}\rho_\mathrm{e}(\boldsymbol{g}_\parallel,g_z)e^{ig_zz}\)と展開します。Hartreeポテンシャルは

\[\begin{equation} V_\mathrm{H}(\boldsymbol{g}_\parallel,z)=\int_{-z_0}^{z_0}dz'\,G(g_\parallel,z,z')\rho_\mathrm{e}(\boldsymbol{g}_\parallel,z')=\sum_{g_z}\rho_\mathrm{e}(\boldsymbol{g}_\parallel,g_z)\int_{-z_0}^{z_0}dz'\,G(g_\parallel,z,z')\,e^{ig_zz'} \label{eq:vh-laue} \end{equation}\]

です。グリーン関数は\(z'\)の指数関数の和なので、積分は解析的に実行できます。bc1（\(g_\parallel\ne0\)）の場合は

\[\begin{equation} V_\mathrm{H}^\mathrm{bc1}(\boldsymbol{g}_\parallel,z)=\sum_{g_z}\frac{4\pi}{g_\parallel^2+g_z^2}\rho_\mathrm{e}(\boldsymbol{g}_\parallel,g_z)e^{ig_zz} -\frac{2\pi}{g_\parallel}\left[e^{g_\parallel(z-z_0)}\sum_{g_z}\frac{\rho_\mathrm{e}(\boldsymbol{g}_\parallel,g_z)e^{ig_zz_0}}{g_\parallel-ig_z}+e^{-g_\parallel(z+z_0)}\sum_{g_z}\frac{\rho_\mathrm{e}(\boldsymbol{g}_\parallel,g_z)e^{-ig_zz_0}}{g_\parallel+ig_z}\right] \label{eq:vh-bc1} \end{equation}\]

となります。第1項は通常の3次元Hartree項\(4\pi\rho_\mathrm{e}(\boldsymbol{G})/G^2\)と同じ形で、第2項は積分範囲が有限であることによる境界項です。bc2とbc3では鏡像電荷項がさらに加わります。たとえばbc3の鏡像電荷項は

\[\begin{equation} V_\mathrm{H}^\mathrm{bc3}(\boldsymbol{g}_\parallel,z)-V_\mathrm{H}^\mathrm{bc1}(\boldsymbol{g}_\parallel,z)=-\frac{2\pi}{g_\parallel}\,e^{g_\parallel(z-2z_1)}\sum_{g_z}\rho_\mathrm{e}(\boldsymbol{g}_\parallel,g_z)\,\frac{2\sinh[(g_\parallel+ig_z)z_0]}{g_\parallel+ig_z} \label{eq:vh-bc3} \end{equation}\]

です。\(g_\parallel=0\)成分は式\(\eqref{eq:g0-bc1}\)–\(\eqref{eq:g0-bc3}\)を用いて同様に求めます。

局所擬ポテンシャル¶

電子–イオン相互作用を、Gauss型の擬イオン電荷\(\rho_\mathrm{g}\)が作る長距離部分、短距離の局所部分、非局所部分に分けます。

\[\begin{equation} E_\mathrm{e\text{-}I}=\iint d\boldsymbol{r}\,d\boldsymbol{r}'\,\rho_\mathrm{e}(\boldsymbol{r})G(\boldsymbol{r},\boldsymbol{r}')\rho_\mathrm{g}(\boldsymbol{r}')+\int d\boldsymbol{r}\,\rho_\mathrm{e}(\boldsymbol{r})V_\mathrm{loc}^\mathrm{short}(\boldsymbol{r})+\sum_i\langle\psi_i|\Delta\hat{V}_\mathrm{PS}|\psi_i\rangle \label{eq:eei} \end{equation}\]

短距離部分と非局所部分は原子核の近くにしか値をもたず、ESMの領域まで届かないので鏡像電荷を生じません。したがって修正は不要です。修正が必要なのは長距離部分\(V_\mathrm{loc}^\mathrm{long}(\boldsymbol{r})=\int d\boldsymbol{r}'\,G(\boldsymbol{r},\boldsymbol{r}')\rho_\mathrm{g}(\boldsymbol{r}')\)だけです。

原子\(\nu\)（価電荷\(Z_\nu\)、位置\(\boldsymbol{\tau}_\nu\)）のGauss電荷を\(\rho_\mathrm{g}^\nu(\boldsymbol{r})=-Z_\nu(\alpha/\pi)^{3/2}e^{-\alpha|\boldsymbol{r}-\boldsymbol{\tau}_\nu|^2}\)とします（電子の電荷を正にとる慣例のため符号は負）。これは面内のGauss関数と\(z\)方向のGauss関数の積なので、\(z'\)積分は解析的に実行でき、\(g_\parallel\ne0\)について

\[\begin{equation} V_\mathrm{loc}^\mathrm{long}(\boldsymbol{g}_\parallel,z)=-\frac{1}{S_0}\sum_\nu Z_\nu\,e^{-i\boldsymbol{g}_\parallel\cdot\boldsymbol{\tau}_{\nu\parallel}}\left[F_\alpha(g_\parallel,z-\tau_{\nu z})+G^\mathrm{i}(g_\parallel,z,\tau_{\nu z})\right] \label{eq:vloc-long} \end{equation}\]

\[\begin{equation} F_\alpha(g_\parallel,\zeta)=\frac{\pi}{g_\parallel}\left[e^{g_\parallel\zeta}\,\mathrm{erfc}\!\left(\frac{g_\parallel}{2\sqrt{\alpha}}+\sqrt{\alpha}\,\zeta\right)+e^{-g_\parallel\zeta}\,\mathrm{erfc}\!\left(\frac{g_\parallel}{2\sqrt{\alpha}}-\sqrt{\alpha}\,\zeta\right)\right] \label{eq:f-alpha} \end{equation}\]

が得られます。\(S_0\)は面内単位胞の面積です。\(F_\alpha\)は裸のCoulomb項\(G^\mathrm{bc1}\)をGauss電荷で平均したもの、\(G^\mathrm{i}\)はグリーン関数の鏡像電荷項（bc1では0、bc2では式\(\eqref{eq:g-bc2}\)の第2項、bc3では式\(\eqref{eq:g-bc3}\)の第2項）です。Gauss電荷がESMの領域と重ならない限り、Gauss電荷の鏡像は点電荷の鏡像と一致するので、鏡像電荷項はGauss関数で平均する必要がありません。

短距離部分\(v_\mathrm{loc}^\mathrm{short}=v_\mathrm{loc}-v_\mathrm{loc}^\mathrm{long}\)は、これまでどおり数値的にFourier変換します。ESMを用いる場合は、通常の計算で逆空間に加える長距離部分の解析的なFourier変換を加えず、代わりに式\(\eqref{eq:vloc-long}\)をLaue表示で加えます。

イオン間相互作用（Ewald和）¶

イオンの電荷を点電荷とすると、イオン間相互作用は

\[\begin{equation} E_\mathrm{I\text{-}I}=\sum_{\boldsymbol{R}_\parallel,\boldsymbol{R}'_\parallel}\sum_{\mu\nu}\frac{Z_\mu Z_\nu}{2}\iint d\boldsymbol{r}\,d\boldsymbol{r}'\,\delta(\boldsymbol{r}-\boldsymbol{R}_{\parallel\mu})\,G(\boldsymbol{r}_\parallel-\boldsymbol{r}'_\parallel,z,z')\,\delta(\boldsymbol{r}'-\boldsymbol{R}'_{\parallel\nu})-E_\mathrm{self} \label{eq:eii} \end{equation}\]

です。\(\boldsymbol{R}_{\parallel\mu}=\boldsymbol{R}_\parallel+\boldsymbol{\tau}_\mu\)は単位胞\(\boldsymbol{R}_\parallel\)内の原子\(\mu\)の位置、\(E_\mathrm{self}\)は自己相互作用です。通常のEwald法と同様に、補助的なGauss電荷を足して引くことで、\(E_\mathrm{I\text{-}I}=E_\mathrm{R}+E_\mathrm{G}-E_\mathrm{self}\)と実空間和と逆空間和に分けます。

実空間和\(E_\mathrm{R}\)には「点電荷−Gauss電荷」の差が現れます。両者の鏡像電荷は同じなので鏡像の寄与は打ち消し合い、裸のCoulomb項だけが残ります。つまり実空間和は通常のEwald和と同じです。

\[\begin{equation} E_\mathrm{R}=\frac{1}{2}\sum_{\mu\nu}\sum_{\boldsymbol{R}_\parallel}{}^{'}\,Z_\mu Z_\nu\,\frac{\mathrm{erfc}\left(\sqrt{\alpha}\,|\boldsymbol{R}_\parallel+\boldsymbol{\tau}_\mu-\boldsymbol{\tau}_\nu|\right)}{|\boldsymbol{R}_\parallel+\boldsymbol{\tau}_\mu-\boldsymbol{\tau}_\nu|} \label{eq:ewald-r} \end{equation}\]

和の\('\)は\(\boldsymbol{R}_\parallel+\boldsymbol{\tau}_\mu-\boldsymbol{\tau}_\nu=0\)の項を除くことを表します。逆空間和は面内の逆格子ベクトル\(\boldsymbol{G}_\parallel\)と\(z\)方向の実空間で表され、\(\boldsymbol{G}_\parallel\ne0\)の部分は

\[\begin{equation} E_\mathrm{G}^{(\boldsymbol{G}_\parallel\ne0)}=\frac{1}{2S_0}\sum_{\mu\nu}Z_\mu Z_\nu\sum_{\boldsymbol{G}_\parallel\ne0}e^{i\boldsymbol{G}_\parallel\cdot(\boldsymbol{\tau}_{\mu\parallel}-\boldsymbol{\tau}_{\nu\parallel})}\left[F_\alpha(G_\parallel,\tau_{\mu z}-\tau_{\nu z})+G^\mathrm{i}(G_\parallel,\tau_{\mu z},\tau_{\nu z})\right] \label{eq:ewald-g} \end{equation}\]

です。\(F_\alpha\)と\(G^\mathrm{i}\)は局所擬ポテンシャルの場合と同じ関数です。\(\boldsymbol{G}_\parallel=0\)の部分は、式\(\eqref{eq:g0-bc1}\)–\(\eqref{eq:g0-bc3}\)の\(|z-z'|\)をGauss関数で平均した

\[\begin{equation} |\zeta|\;\longrightarrow\;\zeta\,\mathrm{erf}(\sqrt{\alpha}\,\zeta)+\frac{1}{\sqrt{\pi\alpha}}\,e^{-\alpha\zeta^2},\qquad \zeta=\tau_{\mu z}-\tau_{\nu z} \label{eq:ewald-g0} \end{equation}\]

に置き換えて、\(E_\mathrm{G}^{(\boldsymbol{G}_\parallel=0)}=\frac{1}{2S_0}\sum_{\mu\nu}Z_\mu Z_\nu\,G(0,\tau_{\mu z},\tau_{\nu z})\)として求めます。自己相互作用項は通常と同じく\(E_\mathrm{self}=\sum_\mu Z_\mu^2\sqrt{\alpha/\pi}\)です。

原子に働く力¶

原子に働く力は全エネルギーの解析的な微分\(\boldsymbol{F}_\mu=-\partial E/\partial\boldsymbol{\tau}_\mu\)として求めます。ESMで変わるのは、グリーン関数を含むEwald項と長距離局所擬ポテンシャル項の寄与だけです。力の計算時間は繰り返しスラブ計算とほぼ同じで、全エネルギーを精度よく保存する分子動力学計算に使えます。

過剰電荷とバイアスの与え方¶

ESM法では、系の電子数を中性の値から変えて過剰電荷を与えます。金属電極を含む境界条件（bc2、bc3）では、過剰電荷を打ち消す対電荷が電極表面に鏡像電荷として誘起され、媒質を含む系全体は自動的に中性に保たれます。一様背景電荷は不要で、静電ポテンシャルの基準（電極の電位＝0）は保たれます。

過剰電荷に対応するバイアスは、帯電した系と中性の系のFermi準位の差として得られます。

\[\begin{equation} \Delta\mu=\mu(\Delta N)-\mu(0) \label{eq:bias} \end{equation}\]

ここで\(\Delta N\)は過剰電子数です。Otani & Sugino（2006）のAl/Si(111)の計算では、\(\Delta\mu\)は誘起電荷にほぼ比例しました。bc2では電荷がスラブの両面に誘起されるため、傾きはbc3のおよそ半分になります。電極から約10 Åの位置にある表面では、1 Vのバイアスでおよそ100表面サイトあたり1個の電子が誘起されました。

bc2では、2つの電極の間に電位差\(V_0\)を直接与えることもできます。境界条件を\(V(\boldsymbol{g}_\parallel,z_1)=0\)、\(V(\boldsymbol{g}_\parallel,-z_1)=V_0\,\delta_{\boldsymbol{g}_\parallel,0}\)とすると、解は式\(\eqref{eq:g-bc2}\)による解\(V^\mathrm{bc2}\)に一様電場の項を加えたものになります。

\[\begin{equation} V(\boldsymbol{g}_\parallel,z)=V^\mathrm{bc2}(\boldsymbol{g}_\parallel,z)-\frac{V_0}{2z_1}(z-z_1)\,\delta_{\boldsymbol{g}_\parallel,0} \label{eq:bias-bc2} \end{equation}\]

電子数を固定してバイアスを結果として得る計算（定電荷計算）は、自己無撞着計算の収束が速く安定です。Otaniら（2008）は、bc3型の配置で過剰電子を与え、バイアスを印加したPt/水界面の第一原理分子動力学計算を行っています。目標とする電位の下で電子数を変化させる計算はConstant-μ法で扱います。

計算コストと適用範囲¶

ESM法の特徴は次のとおりです。

表面・界面に様々な境界条件を課せます。
極性表面や双極子をもつ表面も、補正なしに扱えます。
既存の平面波・擬ポテンシャルコードに容易に導入できます。修正はHartreeポテンシャル、局所擬ポテンシャル、イオン間相互作用の3か所です。
計算時間は従来の方法と同程度です。通常のHartree項の計算は3次元の逆格子ベクトルについて\(O(N^3)\)の演算ですが、ESMでは面内の逆格子ベクトル（\(N^2\)個）ごとに\(g_z\)の和、1次元FFT、\(z\)の和を行うので、演算量は\(N^2(N+N\log N+N)=(2+\log N)N^3\)です。オーダーは変わらず、前因子だけが増えます。

遮蔽媒質は理想金属や真空に限りません。Otani & Sugino（2006）では有限の比誘電率をもつ溶媒と修正Poisson-Boltzmann理論による電解質イオンの分布を組み合わせた境界条件も示されています。また、比誘電率が段差状に変わる元のモデルを滑らかな誘電率分布とDebye-Hückel型の遮蔽に一般化したモデルがHamadaら（2013）によって提案されています。電解液を液体の積分方程式理論で扱う方法はESM-RISM法で扱います。

参考文献¶

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
M. Otani, I. Hamada, O. Sugino, Y. Morikawa, Y. Okamoto, and T. Ikeshoji, "Structure of the water/platinum interface—a first principles simulation under bias potential," Phys. Chem. Chem. Phys. 10, 3609 (2008). 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
I. Hamada, O. Sugino, N. Bonnet, and M. Otani, "Improved modeling of electrified interfaces using the effective screening medium method," Phys. Rev. B 88, 155427 (2013). doi