Constant-μ法¶
Constant-μ法は、電極電位（電子の化学ポテンシャル\(\mu_\mathrm{e}\)）を一定に保って第一原理計算を行う方法です。系の電子数を力学変数とみなし、外部のポテンシオスタットと電子をやり取りさせる仮想的な運動方程式を導入します。各時刻の電子状態計算は電子数を固定した通常の自己無撞着計算のまま行うので、計算の安定性が損なわれません。電極の対電荷はESM法の金属電極が担います。以下では原子単位（\(e=1\)）を用います。
定電荷計算と定電位計算¶
通常の第一原理計算は、系の電子数\(N\)を固定して行います（定電荷計算）。一方、電気化学の実験ではポテンシオスタットで電極電位を制御します（定電位条件）。両者の違いは、電子移動を伴う反応で顕著になります。
例としてPt電極上の水素吸着反応\(\mathrm{H^+}+e^-\to\mathrm{H_{ad}}\)を考えます。電極に過剰電荷を与えた定電荷の分子動力学計算では、電子が電極から吸着種へ移るにつれて、Fermi準位（電極電位）が大きく変化します。実験の条件では、電子はポテンシオスタットから補われ、電極電位は一定に保たれます。定電荷計算は、反応のたびに電極電位が変わってしまう点で実験と対応しません。
グランドポテンシャルによる定式化¶
電子数\(N\)を変数とし、外部から与えた電子の化学ポテンシャル\(\mu_\mathrm{ext}\)の下でのグランドポテンシャル
を考えます。\(E(\{\boldsymbol{r}\};N)\)は原子配置\(\{\boldsymbol{r}\}\)と電子数\(N\)に対するDFTの全エネルギーです。停留条件は
で、系のFermi準位がポテンシオスタットの化学ポテンシャルに一致することを表します。\(\partial^2E/\partial N^2=\partial\mu_\mathrm{e}/\partial N\)は系の微分容量\(C\)の逆数で正なので、この停留点は最小です。電極電位は\(\Phi=-\mu_\mathrm{e}/e\)で定義します。
\(\Omega\)は\(E(N)\)の\(N\)に関するLegendre変換です。定電荷計算だけで\(\Omega\)を求めるには、いくつかの\(N\)で計算し、結果を補間して手作業でLegendre変換する必要があります。
SCFループ内での電子数調整とその問題¶
Lozovoiら（2001）は、自己無撞着場（SCF）計算の中で電子数を調整する方法を提案しました。通常のSCF計算では、各反復でKohn-Sham準位を下から順に電子数が\(N\)になるまで占有します。Lozovoiらの方法では、Fermi準位が\(\mu\)になるまで占有するため、反復ごとに電子数が変わります。
この方法には次の問題があります。
- SCF反復に対して変分原理が成り立たないため、収束が保証されません。通常のSCF計算では全エネルギーが単調に収束しますが、この方法では電荷密度のミキシングをどれだけ小さくしても発散する場合が多くあります。
- 初期電荷密度を注意深く選ぶ必要があり、分子動力学計算のように多数のSCF計算を繰り返す用途には向きません。
したがってSCF計算そのものは定電荷で行い、電子数はSCFループの外側で、Fermi準位を目標値に近づけるように調整する必要があります。そのために、分子動力学の拡張系の方法を用います。
定圧力法との類推¶
Constant-μ法は、Andersen（1980）の定圧力分子動力学法と同じ構成をとります。定圧力法では、セルの体積\(V\)を力学変数とし、原子座標をスケールした座標\(\tilde{\boldsymbol{r}}_i=V^{-1/3}\boldsymbol{r}_i\)で表します。仮想的なラグランジアン
から、体積の運動方程式
が得られます。\(W\)はセルの仮想質量、\(P_\mathrm{in}\)はビリアル定理で与えられる瞬間的な内圧です。内圧が外圧より高ければ体積が増え、低ければ減るというフィードバックが働き、体積は平衡値のまわりを振動して\(\langle P_\mathrm{in}\rangle=P_\mathrm{ext}\)となります。
Constant-μ法では、体積の代わりに電子数を、圧力の代わりに化学ポテンシャルを用います。
| 定圧力法 | Constant-μ法 | |
|---|---|---|
| 示量変数 | 体積\(V\) | 電子数\(N\) |
| 示強変数 | 圧力\(P_\mathrm{in}\) | 化学ポテンシャル\(\mu_\mathrm{in}\) |
| 外部から与える値 | \(P_\mathrm{ext}=\langle P_\mathrm{in}\rangle\) | \(\mu_\mathrm{ext}=\langle\mu_\mathrm{in}\rangle\) |
| 仮想質量 | \(W\) | \(M\) |
| 平衡近傍の復元力 | 体積弾性率\(B=-V_0\,\partial P/\partial V\) | 容量の逆数\(1/C=\partial\mu/\partial N\) |
これまで誰もこの類推を電子数に適用しなかったのは、周期境界条件の下でセルの電荷量を物理的に意味のある形で変える方法がなかったためです（後述）。
仮想電荷粒子の運動方程式¶
ラグランジアンと運動方程式¶
系をポテンシオスタットに接続し、電子数\(N\)を力学変数とします。電子を1個系に移すとポテンシオスタットのエネルギーは\(\mu_\mathrm{ext}\)だけ下がるので、拡張系のポテンシャルエネルギーは\(\tilde{E}=E(\{\boldsymbol{r}\};N)-\mu_\mathrm{ext}N\)です。仮想的なラグランジアンを
とします。\(M\)は電荷量の仮想質量です。Euler-Lagrange方程式は
です。\(\mu_\mathrm{in}=\partial E/\partial N\)は瞬間的な電子系の化学ポテンシャル（Fermi準位）です。Fermi準位が目標より高ければ電子が系から抜け、低ければ電子が入ります。電荷量を表すこの仮想的な自由度を、仮想電荷粒子（fictitious charge particle: FCP）と呼びます。
ハミルトニアンと保存量¶
共役運動量\(\boldsymbol{p}_i=m_i\dot{\boldsymbol{r}}_i\)、\(p_N=M\dot{N}\)を用いると、ハミルトニアンは
となり、正準方程式は
です。\(H_\mu\)は保存量です。Al単原子層（セルあたり4原子、原子位置固定）でFCPの時間発展だけを行ったテスト計算では、Fermi準位は目標値のまわりを単振動し、熱浴を付けない場合の振幅は初期のFermi準位と目標値の差を反映しました。DFTの全エネルギー\(E\)はFermi準位の振動に伴って振動する一方、\(H_\mu\)は保存されました。
Bonnetらの表式と熱浴¶
Bonnetら（2012）は同じ運動方程式を電極電位で表しています。電極電位を\(\Phi=-\partial E/\partial N=-\mu_\mathrm{in}\)、ポテンシオスタットの電位を\(\Phi_\mathrm{pot}=-\mu_\mathrm{ext}\)とすると、FCPに働く力は\(F_N=-\partial\tilde{E}/\partial N=\Phi-\Phi_\mathrm{pot}\)となり、
です。この力は現実の電気的な力と同じ役割を果たし、\(\Phi>\Phi_\mathrm{pot}\)のとき電子を系に引き込みます。平衡化後、\(\Phi\)は\(\Phi_\mathrm{pot}\)のまわりを揺らぎ、\(\langle\Phi\rangle=\Phi_\mathrm{pot}\)となります。
粒子の熱浴（Nosé-Hoover熱浴）だけを付けた場合（方式I）でも、エルゴード性を仮定すれば、温度\(T\)、電位\(\Phi_\mathrm{pot}\)のグランドカノニカル分布が得られます。ただし粒子と電荷の自由度の結合は弱く、平衡化に長い時間がかかります。そこで電荷にも直接熱浴を付けます。粒子と共通の熱浴を使う方式IIと、電荷専用の熱浴（温度\(T_\mathrm{c}\)）を使う方式IIIがあります。方式IIIの運動方程式は、熱浴変数\(\xi_\mathrm{c}\)、その運動量\(P_\mathrm{c}\)と仮想質量\(M_\mathrm{c}\)を用いて
です。\(T_\mathrm{c}=T\)ならグランドカノニカル分布が得られます。\(T_\mathrm{c}<T\)とすると標準的な統計集団からは外れますが、粒子と電荷の結合が弱いため、それぞれ温度\(T\)と\(T_\mathrm{c}\)の準平衡状態が実現し、電極電位の揺らぎを小さくできます。
統計力学的性質¶
熱浴を接続した拡張系は、電子系についてのグランドカノニカル集団を実現します。分配関数は
です。電子数\(N\)の正準分配関数を\(Z_N\)とすると\(\Xi_\mu=\int dN\,Z_N\,e^{\mu_\mathrm{ext}N/k_\mathrm{B}T}\)なので、平均電子数を\(\mu_\mathrm{ext}\)で微分すると、微分容量が電子数の揺らぎで表されます。
これは定圧力法で体積の揺らぎから圧縮率が求まることに対応し、電子回路ではJohnson-Nyquist雑音として知られる関係です。電極電位の標準偏差は\(\sigma_\Phi=\sigma_N/C\propto\sqrt{T/C}\)で、容量は系の大きさに比例するので、小さな系では電位の揺らぎが大きくなります。理想的なポテンシオスタットにつながった電極の電位は揺らがないので、この揺らぎはポテンシオスタットを陽に扱えないことによる人為的なものです。
平衡点のまわりで式\(\eqref{eq:eom-n}\)を線形化すると、\(\mu_\mathrm{in}-\mu_\mathrm{ext}\simeq\delta N/C\)より
となり、電荷量は角振動数\(\omega\)で平均値のまわりを振動します。仮想質量\(M\)は、この振動数が原子の運動と適切に結合するように選びます。Bonnetら（2012）のAl電極の計算では、電荷が\(300\ \mathrm{cm^{-1}}\)（Alのフォノンスペクトルの上端付近）で振動するように\(M\)を選んでいます。
ESM法との組み合わせと電極電位の基準¶
Constant-μ法自体は一般的な方法なので、通常の周期境界条件の下でも実行できます。しかし周期境界条件の下でセルの電荷量を変えると、静電エネルギーの発散を防ぐために一様背景電荷を導入する必要があります。一様背景電荷は便宜的なもので、これを入れるとFermi準位の絶対値が意味を失います。半導体の欠陥の帯電状態の計算でも、絶対値に意味をもたせるには補正が必要です。このため、周期境界条件の下でconstant-μ計算を行っても物理的な意味づけが困難です。
ESM法では、グリーン関数でPoisson方程式を解くため系全体の電荷中性が自然に保たれ、ESM電極に誘起される補償電荷は物理的に意味のある形で分布します。Fermi準位の絶対値の意味も損なわれないので、constant-μ計算の意味づけが明確になります。電気化学セルのモデルとしては、作用極となるスラブの片側にESM金属を対極として置くbc3（真空／スラブ／金属）を用います。
電極電位の基準は次のように決まります。
- bc3ではESM金属の静電ポテンシャルをゼロにとるので、Fermi準位\(\mu_\mathrm{e}\)は対極の電位を基準とした値になります。電極電位\(\Phi=-\mu_\mathrm{e}/e\)は、対極を0 Vとしたときの作用極の電位です。
- 過剰電荷がゼロ（ゼロ電荷電位、pzc）のとき、スラブと対極の間には電場がなく、その間の真空のポテンシャルは対極の電位（ゼロ）に等しくなります。したがって\(\mu_\mathrm{e}\)は対極側表面の仕事関数\(W\)を用いて\(-W\)となり、\(\Phi_\mathrm{pzc}=W/e\)です。
- 電極の相対電位は\(U=\Phi_\mathrm{pot}-\Phi_\mathrm{pzc}\)で定義します。Bonnetら（2012）のAl電極（36原子、対極までの距離10 bohr）では\(\Phi_\mathrm{pzc}=4.31\) Vで、\(U\)を0 V、\(-0.5\) V、\(+0.5\) Vと切り替えると、電極電位は速やかに新しい目標値のまわりの揺らぎに移りました。
電子数\(N\)を連続変数として扱うため、エネルギー関数\(E(\{\boldsymbol{r}\};N)\)が非整数の電子数に対して微分可能である必要があります。近似的な交換相関汎関数では自己相互作用のため、軌道エネルギーが占有数に依存し、非整数の電子数に対する振る舞いが悪くなることがある点に注意が必要です。
構造最適化と反応経路探索¶
分子動力学計算では、各時間ステップで次の手順を繰り返します。
- 現在の原子配置\(\{\boldsymbol{r}\}\)と電子数\(N\)で、定電荷のSCF計算を行う。
- 原子に働く力\(\boldsymbol{F}_i\)とFermi準位\(\mu_\mathrm{in}\)を求める。
- 式\(\eqref{eq:canonical}\)（または式\(\eqref{eq:thermostat}\)）で、原子と電子数を同じ時間刻みで更新する。
電荷の時間発展と電子状態計算が分離しているので、既存の第一原理分子動力学コードにわずかな修正で組み込めます。
同じ考え方で、定電位の下での構造最適化もできます。拡張系のポテンシャルエネルギー\(\tilde{E}=E(\{\boldsymbol{r}\};N)-\mu_\mathrm{ext}N\)に対する一般化された力
を用いて、\(\tilde{E}\)空間で最小化します。電子状態を1回計算するごとに原子とFCP（電荷量）を同時に動かすので、余分なループがなく、通常の構造最適化と同程度のコストで定電位の構造が得られます。
FCPの更新には注意が必要です。最急降下法は最後の収束が遅くなります。2次関数近似による直線探索は、原子を固定した計算ではよく収束しますが、原子を同時に動かすと連続する2ステップで\(F_\mathrm{FCP}\)がほぼ等しくなることがあり、2次関数による次の電荷量の予測が破綻します。電荷量の変化に一定の上限を設ける方法は、セルの表面積が変わると対応できません。そこで、真空のキャパシターのモデルから上限を決めます。電極表面とESMの距離を\(d\)、面積を\(S\)とすると、真空キャパシターの容量は\(C_0=S/4\pi d\)、系の容量は\(C=\varepsilon_\mathrm{sys}S/4\pi d\)（\(\varepsilon_\mathrm{sys}\ge1\)）です。現在の電位差を埋めるのに必要な電荷量の変化は\(C|\mu_i-\mu_\mathrm{ext}|\)ですが、\(C\)は事前にはわかりません。そこで、その下限
を1ステップあたりの電荷量の変化の上限とします。この方法は、金属表面上の分子や電極／固体電解質界面でうまく機能しています。
反応経路の探索には、nudged elastic band（NEB）法を\(\tilde{E}\)空間に適用します（constant-μ+NEB法）。定電荷計算で特定の電位における全エネルギーを得るには、複数の電子数で計算してLegendre変換する必要がありますが、constant-μ法なら望みの電位での値を1回の計算で得られます。反応経路が電位に依存しなければ、定電荷のNEB計算と手作業のLegendre変換で代用できます。しかし現実の系では反応経路が電位に依存するので、constant-μ法が本質的に必要です。Al(001)表面上の水素原子拡散の計算では、グランドポテンシャル\(\Delta\Omega=\Delta E_\mathrm{tot}-\mu_\mathrm{ext}\Delta N\)に沿った経路が得られ、定電荷計算では経路に沿ってFermi準位が変化するのに対し、constant-μ計算では電位を保ったまま経路に沿って電荷量が変化します。
参考文献¶
- 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
- 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
- A. Y. Lozovoi, A. Alavi, J. Kohanoff, and R. M. Lynden-Bell, "Ab initio simulation of charged slabs at constant chemical potential," J. Chem. Phys. 115, 1661 (2001). doi
- H. C. Andersen, "Molecular dynamics simulations at constant pressure and/or temperature," J. Chem. Phys. 72, 2384 (1980). doi