ESM チュートリアル¶

Quantum ESPRESSOに付属するESM法の例題（PW/examples/ESM_example）を実際に計算し、出力される面内平均の電子密度とポテンシャルを図にしながら、ESM法で何が変わるのかを見ていきます。理論の詳細はESM法を参照してください。

準備¶

ESM法は、&system 名前リストの次の2つのキーワードで有効になります。

 &system
    ...
    assume_isolated = 'esm'
    esm_bc = 'bc1'      ! 'pbc', 'bc1', 'bc2', 'bc3'
 /

esm_bc は、スラブ（または分子）の両外側に置く媒質を選びます。

`esm_bc`	左側（\(z<0\)）	右側（\(z>0\)）	主な用途
`pbc`	周期境界条件（ESMを使わない）	同左	比較用
`bc1`	真空	真空	孤立したスラブ・分子
`bc2`	金属	金属	電場をかけた系（コンデンサー）
`bc3`	真空	金属	帯電した電極表面（対極つき）

例題は次のように実行できます。計算はいずれも数秒で終わります。

cd q-e/PW/examples/ESM_example
./run_example_ESM

ここで示す結果は、例題に付属する参照出力（reference/）と全エネルギー・フェルミエネルギーが全桁一致することを確認しています。H₂Oの bc2 と bc3 は例題にない追加計算で、H2O_bc1.in の esm_bc だけを変えたものです。

出力ファイル `*.esm1`¶

ESM法（esm_bc='pbc' を含む）で計算すると、面内（\(xy\)）で平均した電子密度とポテンシャルが prefix.esm1 に書き出されます。

#z (A)  Tot chg (e/A)  Avg v_hartree (eV)  Avg v_local (eV)  Avg v_hart+v_loc (eV)
 -6.00           0.0000048        -198.6208627       198.6194404        -0.0014223

列	内容
1	\(z\)（Å）
2	面内平均した電子密度（e/Å、電子を正にとる。積分すると電子数になる）
3	Hartreeポテンシャル（eV）
4	局所擬ポテンシャル（eV）
5	3列目と4列目の和（eV）

以下の図では5列目、つまり電子が感じる静電的なポテンシャルエネルギーをプロットします。

1. 双極子をもつ分子：H₂O¶

分子の双極子が \(z\) 方向を向くように置いたH₂O（セル 10.58×10.58×12.70 Å、\(\Gamma\) 点）を、4つの境界条件で計算します。

H2O分子の面内平均ポテンシャル

PBC：真空領域でポテンシャルが傾いています。セルの両端で値がそろうという周期境界条件の制約のため、周期的に並んだ双極子がセル全体に人工的な電場をつくります。
bc1（真空/真空）：左右の真空領域がどちらも平坦になり、両側の差は0.424 eVです。これが分子の双極子による真空準位の段差で、孤立した分子の正しい描像です。
bc2（金属/金属）：両端の金属でポテンシャルが0に固定されます。金属に誘起される鏡像電荷が双極子を打ち消すため、真空領域に傾きが生じます。
bc3（真空/金属）：右端の金属で0に固定され、左側の真空は平坦です。

全エネルギーもわずかに変わります。PBCとの差は、周期的に並んだ双極子どうしの相互作用（PBC）や鏡像電荷との相互作用（bc2）の分です。

`esm_bc`	全エネルギー（Ry）	PBCとの差（meV）
`pbc`	−34.43089505	0
`bc1`	−34.43056646	+4.47
`bc2`	−34.43090315	−0.11
`bc3`	−34.43057080	+4.41

2. 金属スラブのフェルミエネルギー：Al(001)¶

Al(001)の1原子層スラブ（4原子、\(6\times6\times1\) の \(k\) 点）の中性状態を計算すると、フェルミエネルギーは次のようになります。

`esm_bc`	フェルミエネルギー（eV）
`pbc`	−2.0409
`bc1`	−4.2462
`bc2`	−4.2503
`bc3`	−4.2474

周期境界条件ではセル内のポテンシャルの平均を0とするため、エネルギーの基準に物理的な意味がありません。ESM法では無限遠の真空（または金属）でポテンシャルが0になるので、フェルミエネルギーがそのまま真空準位を基準にした値、つまり仕事関数の符号を変えたもの（約4.25 eV）になります。ここでの値は、1原子層・低いカットオフの例題用の模型での値です。

3. 帯電した電極：bc3 と `tot_charge`¶

bc3 で系に余分な電荷を与えます（tot_charge = -0.02：電子を0.02個追加、+0.02：0.02個除去）。打ち消す電荷はESMの金属側が担うので、セル全体を一様な背景電荷で中和する必要はありません。

帯電させたAl(001)スラブの電子密度の変化とポテンシャル

電子密度の変化（上）：余分な電荷は、金属（右側）に向いた表面に集まります。左側の真空に向いた表面はほとんど変化しません。
ポテンシャル（下）：スラブと金属の間にだけ電場が生じ（右側の傾き）、左側の真空は平坦なまま全体が上下にずれます。対極（金属）と向かい合った電極表面、すなわちコンデンサーの片側を表しています。

フェルミエネルギーは、中性（−4.25 eV）に対して −0.02e で −3.80 eV、+0.02e で −4.71 eV と、平均0.458 Vずつ変わります。面積（32.8 Å²）あたりの容量に直すと約2.1 μF/cm²で、真空を挟んだ平行平板コンデンサー \(C=\varepsilon_0/d\) に当てはめると \(d\approx4.2\) Åに相当します。これは、スラブの表面から金属までのおおよその距離です。

4. 電場をかける：bc2 と `esm_efield`¶

bc2 で esm_efield = 0.00192148511256006（Ry/bohr、約0.049 V/Å）を与えると、両側の金属の間に一様な電場がかかります。

電場をかけたAl(001)スラブの誘起電荷とポテンシャル

誘起電荷（上）：電子が右側の表面に寄り、左側の表面から減ります。スラブが分極して、金属の内部から電場を遮蔽しています。
ポテンシャル（下）：両側の真空領域に、かけた電場による傾きが現れます。左右の金属の間の電位差は、電場とセルの長さの積（0.049 V/Å × 12 Å ≈ 0.59 V）に対応します。

5. 電極電位を一定に保つ：FCPとGC-SCF¶

ここまでは電荷（tot_charge）を指定して計算しました。実際の電気化学では、電極の電位を指定して、電荷のほうが決まります。これを行うのがConstant-μ法で、目標のフェルミエネルギー \(\mu\) を与えると、系の電荷を自動で調整します。例題（run_example_FCP、run_example_GCSCF）では、中性のときのフェルミエネルギー（零電荷電位、pzc）から ±0.5 Vずれた電位を目標にしています。電極電位 \(U\) とフェルミエネルギーは \(U-U_\mathrm{pzc}=-(\mu-\mu_\mathrm{pzc})/e\) の関係にあるので、+0.5 Vの目標は \(\mu=\mu_\mathrm{pzc}-0.5\) eVです。

Quantum ESPRESSOには2つの方法があります。

! FCP: the charge is relaxed in an outer loop around the SCF
 &control
    calculation = 'relax'
    lfcp = .TRUE.
 /
 &fcp
    fcp_mu = -4.7490          ! target Fermi energy (eV)
    fcp_dynamics = 'bfgs'
    fcp_conv_thr = 1.d-2      ! eV
    freeze_all_atoms = .TRUE. ! move only the charge
 /

! GC-SCF: the charge is adjusted inside the SCF loop
 &control
    calculation = 'scf'
 /
 &system
    lgcscf = .TRUE.
    gcscf_mu = -4.7490        ! target Fermi energy (eV)
    gcscf_conv_thr = 1.d-2    ! eV
 /

FCP（fictitious charge particle）：電荷を仮想的な粒子の座標とみなし、目標とのフェルミエネルギーの差を「力」として、SCFの外側で電荷を動かします。電荷を変えるたびにSCFを解き直します。
GC-SCF（grand-canonical SCF）：SCFの反復の中で、電子密度と一緒に電荷も更新します。1回のSCFで電荷とフェルミエネルギーがそろいます。

bc3 で +0.5 Vと −0.5 Vを目標にしたときの、フェルミエネルギーと電荷の変化を、SCF反復の累計に対して比べます。

FCPとGC-SCFが目標のフェルミエネルギーに近づく様子

FCP（青）は、SCFを1回解き終えるごとに電荷を更新するので、点と点の間に10〜20回程度のSCF反復が入ります。GC-SCF（橙）は、1回のSCFの中で電荷が目標の値へ向かい、10回ほどの反復で収束します。最終的な電荷は両者でほぼ同じです。

境界条件	目標	FCP：SCF回数	FCP：反復の累計	FCP：電荷 (e)	GC-SCF：反復	GC-SCF：電荷 (e)
`bc3`	+0.5 V	5	68	+0.0217	11	+0.0218
`bc3`	−0.5 V	5	63	−0.0220	12	−0.0219
`bc2`	+0.5 V	3	25	+0.0437	10	+0.0438
`bc2`	−0.5 V	5	57	−0.0437	12	−0.0436

この結果から、次のことが読み取れます。

bc3 で0.5 Vあたり約0.0218 eの電荷が出入りします。これは3節で電荷から見積もった容量（0.02 eで0.458 V）と一致します。
bc2 ではその約2倍の電荷が出入りします。スラブの両面がそれぞれ金属と向かい合い、コンデンサーが2つ並列になっているためです。
2つの方法では、反復のかかり方が違います。GC-SCFは1回のSCFで目標の電位まで持っていくので効率的ですが、電子密度と電荷を同時に収束させるため、1回のSCFの反復回数は多めになりがちです（例題でも電荷密度の混合の割合を mixing_beta = 0.1 と小さくしています。この例題は系が小さいため、中性のSCF（12回）と同程度の11〜12回で済んでいます）。一方、FCPの1回のSCFは電子数を固定した通常のSCFなので、1回あたりの反復は短い傾向にあります。ただし、そのSCFを何度も繰り返すため、反復の合計は多くなります。

これらの特徴から、構造を固定してフェルミエネルギーを目標に合わせる計算にはGC-SCFが有利です。一方、FCPは電荷を運動方程式で扱うので、分子動力学と組み合わせて、電位を一定に保ったまま原子を動かす計算に有利です。分子動力学では前のステップの電荷と電子密度がよい初期値になるため、各ステップの電子数一定のSCFは短く済みます（理論はConstant-μ法）。

まとめ¶

知りたいこと	境界条件	ESMで得られること
双極子をもつ系・非対称なスラブ	`bc1`	人工的な電場のない孤立系と、真空準位の段差
仕事関数・エネルギー準位の絶対値	`bc1`〜`bc3`	真空（または金属）を基準にしたフェルミエネルギー
帯電した電極表面	`bc3` ＋ `tot_charge`	対極を含んだ電荷の分布と容量
電場をかけた系	`bc2` ＋ `esm_efield`	分極と遮蔽
電位を指定した電極	`bc2`・`bc3` ＋ FCP／GC-SCF	電位に応じた電荷と容量