ESM-RISM チュートリアル¶

参照相互作用サイトモデル（RISM）を組み合わせたDFT計算のチュートリアル

このチュートリアルでは、QUANTUM ESPRESSO に実装された、密度汎関数理論（DFT）と参照相互作用サイトモデル（RISM）を組み合わせた計算の使い方を学びます。

DFT計算では、電子状態を量子力学に基づいて扱います。一方、RISMは古典的な溶液理論で、溶液を古典力場で扱います。
したがって、DFT+RISMを用いると、溶液中の原子、イオン、分子、物質の電子状態を計算できます。この枠組みでは、溶媒の電荷密度は、分子内構造を固定した陰溶媒として扱われます。
さらに、DFT+RISMの大きな利点は、溶液中のイオンや分子に対するグランドカノニカル条件のもとで、電極/電解液界面を精密に扱えることです。

理論的背景はS. Nishihara and M. Otani, Phys. Rev. B 96, 115429 (2017)に簡潔にまとめられています。理論の考え方を理解しておくとRISM計算を学ぶうえで大いに役立つので、まずこの論文を読んでください。

PWscfには3種類のRISM計算法が実装されています。

このチュートリアルは次の内容で構成されています。

RISM計算の入力値とファイル
NaCl水溶液中のH\(_2\)O分子に対する3D-RISM計算
NaCl水溶液中のAl(111)スラブに対するESM-RISM計算
ESM-RISM計算による標準水素電極（SHE）電位の決定

0. RISM計算の入力値¶

DFT+RISM計算を始める前に、RISM計算の入力値とファイルを示します。
RISM計算を行う際は、次の一覧と説明を参照してください。

RISM計算の入力値
*.MOL ファイル

0-1. RISM計算の入力値¶

ここでは、RISM計算でよく使う入力値を簡単に説明します。
まず、入力値とその説明を読んでください。

`Namelist: &CONTROL`¶

trism: LOGICAL/default = .FALSE.

RISM計算のスイッチです。この値を設定しない場合、RISM関連の入力値は無視されます。この値を有効にした場合は、&RISM 名前リストと SOLVENTS カードが必要です。RISM計算では、溶媒の分布をRISM理論で計算し、陽に扱う溶質（表面スラブなど）をDFTで扱います。
assume_isolated = esm かつ esm_bc=bc1 とした場合は、Laue-RISM 計算が行われます。

mixing_beta のデフォルト値は、3D-RISMとLaue-RISMのいずれでも0.2です。

BFGSによる構造緩和では、ignore_wolfe は常に .TRUE. になります。

`Namelist: &RISM`¶

RISM計算の主な入力値を記述する名前リストです。trism=.TRUE. のときに使用できます。

nsolv: INTEGER/default = 0

単位胞内の溶媒（すなわち分子種）の数です。RISM計算では必ず指定する必要があります。

closure: CHARACTER/default = 'kh'

閉包式の種類を指定します。
現在、次の2種類の閉包式を使用できます。

'kh': Kovalenko–Hirataモデル [A.Kovalenko, F.Hirata, JCP 110, 10095 (1999)]
'hcn': HyperNetted-Chainモデル [J.P.Hansen et al., Theory of simple liquids. Academic Press, London, 1990]

注：'kh' 閉包式を強く推奨します。

tempv: REAL/default = 300.0

溶媒の温度を指定します（単位はケルビン）。

ecutsolv: REAL/default = 4 × ecutwfc

溶媒の相関関数に対する運動エネルギーのカットオフ（単位はRy）です。溶質が孤立系やスラブの場合は、デフォルト値を使ってかまいません。骨格構造をもつ溶質や多孔質の溶質（ゼオライト、MOFなど）では、より大きな値を用いることが望まれます。骨格内に閉じ込められた溶媒は、しばしば高い周波数成分をもつためです。

solute_lj(1:ntyp): CHARACTER/default = 'uff'

ityp 番目の溶質の原子種について、溶質–溶媒相互作用のLennard-Jonesポテンシャルの種類を指定します。

'uff': Universal Force Field [A.K.Rappe et al., JACS 144, 10024 (1992)]
'clayff': Clay's Force Field [R.T.Cygan et al., JPC B 108, 1255 (2004)]
'opls-aa': OPLS-AA（QM/MM用の汎用パラメータ）[https://github.com/thatchristoph/vmd-cvs-github/blob/master/plugins/bossconvert/oplsaa.par を参照]
'none': Lennard-Jonesパラメータを、ユーザーが solute_epsilon と solute_sigma で指定します。

注：この値は3D-RISMまたはLaue-RISM計算で使用されます。

solute_epsilon(1:ntyp): REAL/default =

ityp 番目の原子種のLennard-Jonesパラメータ \(\varepsilon\)（単位はkcal/mol）を指定します。
solute_lj(1:ntyp)='none' の場合は、この値を必ず設定してください。

solute_sigma(1:ntyp): REAL/default =

ityp 番目の原子種のLennard-Jonesパラメータ \(\sigma\)（単位はÅ）を指定します。
solute_lj(1:ntyp)='none' の場合は、この値を必ず設定してください。

starting1d: CHARACTER/default = 'zero'

1D-RISM計算で求める相関関数の初期値を指定します。

'zero': 初期の相関関数を1D-RISMで新たに計算します。
'file': 初期の相関関数を既存の 1d-rism_csvv_r.xml ファイルから読み込みます。
リスタート計算では、これがデフォルトです。
'fix': 初期の相関関数を既存の "1d-rism_csvv_r.xml" ファイルから読み込み、
nscf と bands の計算では1D-RISMを計算しません。

starting3d: CHARACTER/default = 'zero'

3D-RISMまたはLaue-RISM計算で求める相関関数の初期値を指定します。

'zero': 初期の相関関数を3D-RISMで新たに計算します。
'file': 初期の相関関数を既存の 3d-rism_csuv_r.dat ファイルから読み込みます。
リスタート、nscf、bands の計算では、これがデフォルトです。

smear1d: REAL/default = 2.0

1D-RISMのクーロン相互作用のスメアリング半径です。
注：通常、この値は変更しません。

smear3d: REAL/default = 2.0

3D-RISMとLaue-RISMのクーロン相互作用のスメアリング半径です。
注：通常、この値は変更しません。

rism1d_maxstep: INTEGER/default = 50000

1D-RISMの反復の最大回数です。
注：rism1d_maxstep=10000 以上を推奨します。

rism3d_maxstep: INTEGER/default = 5000

3D-RISMとLaue-RISMの反復の最大回数です。
注：rism3d_maxstep=10000 以上を推奨します。

rism1d_conv_thr: REAL/default = 1.d-8

1D-RISMの収束判定条件です。
注：デフォルト値を推奨します。

rism3d_conv_thr: REAL/default = 1.d-5

3D-RISMとLaue-RISMの収束判定条件です。
注：フェルミ準位が収束判定条件に依存することがわかっているため、rism3d_conv_thr=1.d-6 を推奨します。

mdiis1d_size INTEGER/default = 20

1D-RISMのModified DIIS（MDIIS）のサイズです。
MDIISは、収束した相関関数を得るための頑健な最適化手法の一つです。

mdiis3d_size INTEGER/default = 10

3D-RISMとLaue-RISMのMDIISのサイズです。
注：mdiis3d_size \(\geq\) 10 を推奨します。

mdiis1d_step REAL/default = 0.5

1D-RISMのMDIISのステップ幅です。
注：mdiis1d_step \(=\) 0.1 を推奨します。
1D-RISM計算が10000ステップを超えても収束しない場合は、
mdiis1d_step を0.1より小さくするか、0.1から0.12に大きくしてみてください。

mdiis3d_step REAL/default = 0.8

3D-RISMとLaue-RISMのMDIISのステップ幅です。
注：3D-RISMまたはLaue-RISM計算が収束しない場合は、
mdiis3d_step を変えてみてください。

rism1d_bond_width REAL/default =

1D-RISMで分子内相関をスメアするための、結合のガウス幅です。

3D-RISM：デフォルト = 0.0
Laue-RISM：デフォルト = \(2/\sqrt{\rm ecutwfc}\)

rism1d_nproc INTEGER/default = 128

1D-RISMの計算に使うプロセス数です。
FFTグリッドに使うプロセス数が \(\leq 128\) の場合、rism1d_nproc は適切な値に自動で決まります。

rism3d_conv_level REAL/default = 0.1

3D-RISMの収束レベルです。
lgcscf=.TRUE. の場合、rism3d_conv_level のデフォルト値は0.3になります。

rism3d_conv_level は [0.0, 1.0] の範囲で設定する必要があります。

rism3d_conv_level = 0.0：収束レベルは 'low' です。
推定されるエネルギーの誤差が conv_thr より十分大きいとき、3D-RISMまたはLaue-RISMの収束判定条件は rism3d_conv_thr より大きくなります。
0.0 < rism3d_conv_level < 1.0：収束レベルは 'medium' です。
3D-RISMまたはLaue-RISMの収束判定条件は 'low' と 'high' の中間の値で、rism3d_conv_level がその混合比です。
rism3d_conv_level = 1.0：収束レベルは 'high' です。
3D-RISMまたはLaue-RISMの収束しきい値は常に rism3d_conv_thr に等しくなります。

注：3D-RISMやLaue-RISMの反復法の振る舞いは、溶液の系に強く依存します。RISM-SCF計算が収束しないこともあります。その場合、rism3d_conv_level は収束を得るのに非常に有効です。
私たちの経験では、3D-RISMやLaue-RISMには rism3d_conv_level = 0.5または1.0がよい値です。

rism3d_planar_average LOGICAL/default =

面内平均した溶媒と溶質の電荷密度および静電ポテンシャルの出力を制御します。
Laue-RISM計算では rism3d_planar_average は .TRUE. となり、これらの量は outdir/prefix.rism1 ファイルに書き出されます。

laue_nfit INTEGER/default = 4

ESM計算と同様に、セルの端で多項式フィットに使うz方向のグリッド点の数です。Laue-RISMでのみ使用されます。

laue_expand_right REAL/default = -1.0

Laue-RISMでのみ使用されます。laue_expand_right に正の値を与えると、単位胞の右側の溶媒領域の終端位置を、単位胞の端からのオフセット（単位はbohr）として設定します（溶媒領域は z = + [L_z/2 + laue_expand_right] で終わります）。ここで、L_zは単位胞のz方向の長さです。

laue_expand_left REAL/default = -1.0

Laue-RISMでのみ使用されます。laue_expand_left に正の値を与えると、単位胞の左側の溶媒領域の終端位置を、単位胞の端からのオフセット（単位はbohr）として設定します（溶媒領域は z = - [L_z/2 + laue_expand_left] で終わります）。ここで、L_zは単位胞のz方向の長さです。

laue_starting_right REAL/default = 0.0

Laue-RISMでのみ使用されます。単位胞の右側の溶媒領域の開始位置（単位はbohr）を設定します。溶媒領域は [ laue_starting_right , L_z/2 + laue_expand_right ] と定義され、この範囲で分布関数が有限の値をとります。ここで、L_zは単位胞のz方向の長さです。

laue_starting_left REAL/default = 0.0

Laue-RISMでのみ使用されます。単位胞の左側の溶媒領域の開始位置（単位はbohr）を設定します。溶媒領域は [ -(L_z/2 + laue_expand_left), laue_starting_left ] と定義され、この範囲で分布関数が有限の値をとります。ここで、L_zは単位胞のz方向の長さです。

laue_buffer_right REAL/default =

Laue-RISMでのみ使用されます。
デフォルト値は次のように設定されます。

laue_expand_right > 0.0 の場合、laue_buffer_right = 8.0
laue_expand_right \(\leq\) 0.0 の場合、laue_buffer_right = -1.0

正の値の場合、単位胞の右側の溶媒領域のバッファ長（単位はbohr）を設定します。このとき、相関関数は [laue_starting_right-laue_buffer_right, L_z/2 + laue_expand_right ] の内側で定義されます。ここで、L_zは単位胞のz方向の長さです。

laue_buffer_left REAL/default =

Laue-RISMでのみ使用されます。
デフォルト値は次のように設定されます。

laue_expand_left > 0.0 の場合、laue_buffer_left = 8.0
laue_expand_left \(\leq\) 0.0 の場合、laue_buffer_left = -1.0

正の値の場合、単位胞の左側の溶媒領域のバッファ長（単位はbohr）を設定します。このとき、相関関数は [-(L_z/2 + laue_expand_left), laue_starting_left+laue_buffer_left] の内側で定義されます。ここで、L_zは単位胞のz方向の長さです。

`Card: SOLVENTS`¶

RISM計算で、溶媒分子、溶媒のMOLファイル、溶媒密度を指定する省略可能なカードです。
trism が .FALSE. の場合、このカードは無視されます。

書式：

SOLVENTS {Density_unit}
Molecule_name(1), Density(1), MOL_file(1)
....
Molecule_name(nsolv), Density(nsolv), MOL_file(nsolv)

各項目の説明：

Density_unit: CHARACTER

1/cell：溶媒密度を、単位胞内の分子数で定義します。
mol/L：溶媒密度を、モル濃度で定義します。
g/cm^3：溶媒密度を、1 cm\(^3\) あたりのグラム数で定義します。

Molecule_name: CHARACTER

溶媒分子のラベルです。

density: REAL

溶媒分子の密度（単位は Density_unit）です。density が負の場合、溶媒分子の密度は指定したMOLファイルから読み込まれます。

MOL_file: CHARACTER

溶媒分子の .MOL ファイルです。.MOL ファイルについては次節で説明します。

0-2. MOLファイル¶

.MOL ファイルには、分子やイオンの密度、分子やイオンの電荷、xyzファイル形式の分子構造、各原子のLennard-Jonesパラメータなどが含まれます。RISM計算を実行するには、.MOL ファイルを、擬ポテンシャルファイルのパスである pseudo_dir/ に置く必要があります。

このチュートリアルに必要な .MOL ファイルと擬ポテンシャルファイルは、./MOL/ ディレクトリにあります。
新しい溶質や溶媒の分子・イオンを使いたい場合は、Satomichi NishiharaのGitHubリポジトリ [https://github.com/nisihara1/MOLs] から入手してください。必要なMOLファイルが見つからない場合は、目的に合わせて新しいMOLファイルを作成してみてください。

1. NaCl水溶液中のH₂O分子に対する3D-RISM計算¶

この節では、陽に扱う溶質としてのH\(_2\)O分子が、陰溶媒としてのNaCl水溶液（1 mol/L）中にある系について、3D-RISM計算を行います。
このチュートリアルでは次の入力ファイルを使います。まず、入力ファイルと、必要な *.MOL ファイルおよび擬ポテンシャルファイルを計算ディレクトリにコピーしてください（.MOL ファイルと擬ポテンシャルファイルは、それぞれ ./MOL/ と ./pseudo/ ディレクトリにあります）。

この計算では、水分子にSPCモデルのLennard-Jones（LJ）パラメータを使います。
RISM計算は一種のモデル計算なので、入力パラメータの扱いには常に注意が必要です。
計算の前にLJパラメータを確認したい場合は、.MOL ファイルを開いて読んでください。

次に、エディタで入力ファイルを開き、pseudo_dir を擬ポテンシャルファイルのパスに書き換えてから、入力を読んで計算の内容を理解してください。入力の意味を忘れた場合は、上の入力値の一覧をもう一度確認してください。

入力ファイルを読んで理解できれば、3D-RISM計算の準備は完了です。pw.x を使って、最初の3D-RISM計算を実行してみましょう。

&control
   calculation  = 'relax'
   prefix       = 'H2O.NaCl_aq.1mol'
   pseudo_dir   = './pseudo'
   outdir       = './output'
   trism        = .true.
   nstep        = 100
/
&system
   ibrav = 1, a = 20.0
   ntyp = 2, nat = 3
   ecutwfc = 40.0, ecutrho = 320.0
   nbnd = 8
   occupations = 'smearing', smearing = 'gauss', degauss = 0.01
/
&electrons
   electron_maxstep = 200
/
&ions
   ion_dynamics = 'bfgs'
/
&rism
   nsolv    = 3, closure  = 'kh'
   tempv    = 300.0  ! Kelvin
   ecutsolv = 160.0  ! Rydberg
   !
   ! Lennard-Jones for each atom (SPC model)
   !
   ! For O (TIP5P)
   solute_lj(1) = 'none'
   solute_epsilon(1) = 0.1554  ! kcal/mol
   solute_sigma(  1) = 3.1660  ! angstrom
   ! For H (modified H)
   solute_epsilon(2) = 0.0460
   solute_sigma(  2) = 1.0000
   !
   ! 1D-RISM's setting
   !
   starting1d      = 'zero'
   rism1d_conv_thr = 1.0e-8
   rism1d_maxstep  = 10000
   mdiis1d_size    = 20
   mdiis1d_step    = 0.5
   !
   ! 3D-RISM's setting
   !
   starting3d      = 'zero'
   rism3d_maxstep  = 10000
   rism3d_conv_thr = 1.0e-6
/
ATOMIC_SPECIES
 O  -1.0  O.pbe-n-rrkjus_psl.1.0.0.UPF
 H  -1.0  H.pbe-rrkjus_psl.1.0.0.UPF
ATOMIC_POSITIONS {angstrom}
 O  10.00000  10.00000  10.00000  0 0 0
 H  10.81649  10.57736  10.00000
 H   9.18351  10.57736  10.00000

K_POINTS {gamma}

SOLVENTS {mol/L}
 H2O  -1.0  H2O.spc.MOL ! Now we use density of H2O
 Na+   1.0  Na+.aq.MOL
 Cl-   1.0  Cl-.aq.MOL

RISM計算が正しく終了すると、1D-RISMと3D-RISMの計算結果が得られます。

まず、計算のログファイルを開くと、1D-RISMで求めた溶媒和の化学ポテンシャルの結果があります。
各分子・イオンについての合計の値が、溶媒和自由エネルギーに対応します。

次に、ログファイルには3D-RISM計算のログもあります。
ここでも、溶媒和の化学ポテンシャル、溶媒の総数と総電荷、DFT+RISM計算の全エネルギーなどを確認できます。

この計算では、陽に扱う水分子の構造最適化も行っています。原子にはたらく力を確認してください。計算全体が収束に近づくにつれて、力は小さくなります。

次に、1D-RISM計算の結果をさらに詳しく見てみます。計算後、/WORKDIR/output/ ディレクトリに H2O.NaCl_aq.1mol.1drism ファイルができているはずです。
このファイルには、溶液を構成する原子の組ごとの対分布関数のデータが含まれています。
下の図に、すべての対分布関数を示します。これらの分布関数は、\(r\) が十分大きい極限で1に規格化されていることに注意してください。つまり、各分子の密度は掛けられていません。

NaCl水溶液中のH2Oについて1D-RISMで得られた対分布関数

例として、1D-RISMの結果を用いて、着目する原子のまわりの原子数を計算します。
原子数は、密度と分布関数から次のように評価できます。

\[ \begin{aligned} N_{\rm at}= \int 4\pi r^2 \rho g(r) dr \end{aligned} \]

下の図に、動径分布関数と原子数を示します。青と赤の曲線は、それぞれO–Na\(^+\) とNa\(^+\)–Cl\(^-\) の組の分布関数です。橙と緑の曲線は、それぞれNa\(^+\) イオンのまわりのO原子とCl\(^-\) イオンの数を表します。分布関数の第1ピークの位置に、それぞれ1個のO原子と1個のCl\(^-\) イオンが見いだされます。

Na+イオンのまわりの動径分布関数と原子数

さらに解析を進めるために、3D-RISMで得た溶媒の電荷分布も可視化します。
3次元の密度分布ファイルを作るには、実行ファイル pprism.x を使います。
pprism.x には、pp.x と似た形式の入力ファイルが必要です。pprism.x の入力ファイルを以下に示します。

&inputpp
  outdir  = './output'
  prefix  = 'H2O.NaCl_aq.1mol'
  lpunch  = .true.
/
&plot
  output_format = 6
  iflag = 3
  nx = 160, ny = 160, nz = 160
  x0(1) = 0.0
  x0(2) = 0.0
  x0(3) = 0.0
  interpolation = 'fourier'
/

pprism.x を実行すると、次のcubeファイルが得られます（output_format=6）。

H2O.NaCl_aq.1mol.3drism_chg.cube
H2O.NaCl_aq.1mol.3drism_rho#1_O.cube
H2O.NaCl_aq.1mol.3drism_rho#2_H.cube
H2O.NaCl_aq.1mol.3drism_rho#3_Na.cube
H2O.NaCl_aq.1mol.3drism_rho#4_Cl.cube
H2O.NaCl_aq.1mol.3drism_v_solu.cube
H2O.NaCl_aq.1mol.3drism_v_solv.cube
H2O.NaCl_aq.1mol.3drism_v_tot.cube

ここでは、H2O.NaCl_aq.1mol.3drism_chg.cube を VESTA で描画した結果を示します。

陽に扱ったH2O分子のまわりの溶媒の電荷分布の等高線図

陽に扱ったH\(_2\)O分子のまわりの溶媒の電荷分布の等高線図。

H\(_2\)O分子のO原子は負に帯電しているので、溶媒の正電荷（赤）がクーロン相互作用によってO原子に引き寄せられます。
一方、H原子は正に帯電しています。このため、溶媒の負電荷（青）はH原子の近くに分布します。
さらに、溶質から遠い領域では、電荷分布は中性（緑）になります。この領域では、クーロン相互作用が溶媒によって遮蔽されていることがわかります。

以上のように、陽に扱った溶質のまわりの溶媒の電荷分布を可視化すると、溶媒–溶質相互作用を理解するのに非常に役立ちます。

2. NaCl水溶液中のAl(111)スラブに対するESM-RISM計算¶

ここでは、次の5つの計算を行います。

2-1. NaCl水溶液中の電荷中性のAl(111)スラブに対するESM-RISM計算
2-2. NaCl水溶液中の正に帯電したAl(111)スラブに対するESM-RISM計算
2-3. NaCl水溶液中の負に帯電したAl(111)スラブに対するESM-RISM計算
2-4. NaCl水溶液中で溶媒イオン（Na\(^+\)）を陽に扱ったAl(111)スラブに対するESM-RISM計算
2-5. NaCl水溶液中で溶媒イオン（Cl\(^-\)）を陽に扱ったAl(111)スラブに対するESM-RISM計算

これらはいずれもESM-RISM計算です。
ESM-RISMでは、開放境界条件を使うために assume_isolated = "esm" と esm_bc = "bc1" が必要です。
したがって、面内方向には周期境界条件を、表面に垂直な方向には開放境界条件を用います。

ここで、ESM-RISM計算の利点を簡単に説明します。ESM-RISMの枠組みはグランドカノニカル条件を扱えるため、帯電した表面を扱えます。通常、周期境界条件を用いるDFTの枠組みでは、静電ポテンシャルの項が発散するため、帯電した表面を扱えません。発散を避けるために、通常は一様背景電荷（ジェリウム）を導入します。しかし、この場合は全エネルギーの原点を一意に定義できません。このため、帯電した表面と溶液の界面を扱うことは困難でした。

これに対して、ESM-RISMでは、溶質と溶媒からなる系全体で電荷中性条件が常に満たされます。
このため、ESM-RISMは、電気化学界面のような物質/溶液界面を調べるのに非常に有用な手法です。

2-1. NaCl水溶液中の電荷中性のAl(111)スラブに対するESM-RISM計算¶

まず、Al(111)/NaCl水溶液界面を例に、ESM-RISM計算の実行方法を学びます。
入力ファイルと入力値を読んで、ESM-RISM計算の内容を理解しましょう。&system には、bc1（真空/スラブ/真空）のESM計算の記述があります。つまり、z方向には開放境界条件、xy面内には周期境界条件という混合境界条件を用います。ESM-RISM計算でも、ATOMIC_POSITIONS の定義は通常のESM計算と同じであることに注意してください（つまり、座標の原点は単位胞の中心に移されます）。
&rism には、新しい入力値 laue_expand_right があります。この値は、ESM-RISMで単位胞を拡張する長さを決めます。ESM-RISMでは、表面から遠い位置で溶媒の密度分布が一定になる必要があります。そのためには十分な厚さの真空層が必要ですが、厚い真空層をとったDFT計算には膨大な計算コストがかかります。これを避けるために、laue_expand_right（または laue_expand_left）による拡張セルの手法を用います。laue_expand_right の詳細は、RISM計算の入力値の一覧を参照してください。

擬ポテンシャルファイルのパス（pseudo_dir）を環境に合わせて変更し、pw.x を実行してください。

&control
   calculation  = 'scf'
   prefix       = 'Al111.NaCl_aq'
   pseudo_dir   = './pseudo'
   outdir       = './output'
   trism        = .true.
/
&system
   ibrav = 0
   ntyp = 1, nat = 3
   ecutwfc = 25.0, ecutrho = 200.0
   nbnd = 8
   occupations = 'smearing', smearing = 'gauss', degauss = 0.01
   assume_isolated = "esm", esm_bc = "bc1", tot_charge = 0.0
/
&electrons
   electron_maxstep = 200
   mixing_beta = 0.1
   diagonalization = 'rmm', diago_rmm_conv = .false., diago_rmm_ndim = 4
/
&rism
   nsolv    = 3, closure  = 'kh'
   tempv    = 300.0  ! Kelvin
   ecutsolv = 160.0  ! Rydberg
   !
   ! Lennard-Jones for each atom
   !
   ! For Al using universal force field
   solute_lj(1) = 'uff'
   !
   ! 1D-RISM's setting
   !
   starting1d      = 'zero'
   rism1d_conv_thr = 1.0e-8
   rism1d_maxstep  = 10000
   mdiis1d_size    = 20
   mdiis1d_step    = 0.5
   !
   ! 3D-RISM's setting
   !
   starting3d      = 'zero'
   rism3d_maxstep  = 10000
   rism3d_conv_thr = 1.0e-6
   !
   ! For Laue-RISM calculation
   ! in a.u.
   !
   laue_expand_right = 100.0
/
ATOMIC_SPECIES
 Al  -1.0  Al.pbe-nl-rrkjus_psl.1.0.0.UPF
CELL_PARAMETERS angstrom
2.85671139599365 0.00000000000000 0.00000000000000
1.42835569799683 2.47398464021101 0.00000000000000
0.00000000000000 0.00000000000000 35.0000000000000
ATOMIC_POSITIONS angstrom
Al 1.4283556980 0.8246615467 -2.3324950875
Al 0.0000000000 1.6493230935 -0.0000000000
Al 0.0000000000 0.0000000000  2.3324950875
K_POINTS {automatic}
 8 8 1 1 1 0
SOLVENTS {mol/L}
 H2O  -1.0  H2O.spc.MOL
 Na+   1.0  Na+.aq.MOL
 Cl-   1.0  Cl-.aq.MOL

ESM-RISM計算が終わると、output/ ディレクトリに prefix.rism1 ファイルができます。このファイルには、ESM-RISM計算の重要な結果が含まれています。prefix.rism1 ファイルには、面内平均した静電ポテンシャルと溶液の電荷密度の情報が書かれています。以下では、prefix.rism1 ファイルに書かれた結果を描画します。

下の図に、溶液の静電ポテンシャル（左）と、NaCl水溶液の溶媒の電荷密度（中央と右）を示します。

静電ポテンシャル¶

溶媒のポテンシャル（黒線）は、Lennard-Jonesポテンシャルと静電ポテンシャルの和を介した表面–溶媒相互作用のため、Al表面の近くで振動しています。この振動は全ポテンシャル（青線）にも見られます。
開放境界条件では、エネルギーの原点は左右の真空準位で定義されます。しかし、ESM-RISMでは、エネルギーの原点を溶媒のバルク領域に移します。

溶媒密度¶

溶媒密度は、スラブの中心から測って約4.5 Åから \(z\rightarrow \infty\) まで分布しています。さらに、溶媒密度は表面近くから約20 Åまで振動し、その先で一定値に収束します。この振動は、LJポテンシャルと静電ポテンシャルの和を介した表面と溶媒の相互作用によるものです。密度分布が一定であることは、溶媒が一様に分布していることを示しています。
加えて、溶媒密度が一様に分布していることにより、グランドカノニカル条件での計算が可能になります。つまり、一様な密度の領域が溶媒のリザーバーの役割を果たします。これはESM-RISM計算の重要な利点です。以下の帯電した表面の計算で、この利点がわかります。

NaCl水溶液中の電荷中性のAl(111)スラブの静電ポテンシャルと溶媒密度

2-2. NaCl水溶液中の正に帯電したAl(111)スラブに対するESM-RISM計算¶

次に、正に帯電したAl(111)表面とNaCl水溶液の界面を計算します。この計算では、tot_charge を用いて表面電荷 \(\sigma=+0.02e\) を与えます。その他の入力は2-1節と同じです。

&control
   calculation  = 'scf'
   prefix       = 'Al111-0.02e.NaCl_aq'
   pseudo_dir   = './pseudo'
   outdir       = './output'
   trism        = .true.
/
&system
!   ibrav = 1, a = 10.0
   ibrav = 0
   ntyp = 1, nat = 3
   ecutwfc = 25.0, ecutrho = 200.0
   nbnd = 8
   occupations = 'smearing', smearing = 'gauss', degauss = 0.01
   assume_isolated = "esm", esm_bc = "bc1", tot_charge = 0.02
/
&electrons
   electron_maxstep = 200
   mixing_beta = 0.1
   diagonalization = 'rmm', diago_rmm_conv = .false., diago_rmm_ndim = 4
/
&rism
   nsolv    = 3, closure  = 'kh'
   tempv    = 300.0  ! Kelvin
   ecutsolv = 160.0  ! Rydberg
   !
   ! Lennard-Jones for each atom
   !
   ! For Al using universal force field
   solute_lj(1) = 'uff'
   !
   ! 1D-RISM's setting
   !
   starting1d      = 'zero'
   rism1d_conv_thr = 1.0e-8
   rism1d_maxstep  = 10000
   mdiis1d_size    = 20
   mdiis1d_step    = 0.5
   !
   ! 3D-RISM's setting
   !
   starting3d      = 'zero'
   rism3d_maxstep  = 10000
   rism3d_conv_thr = 1.0e-6
   !
   ! For Laue-RISM calculation
   ! in a.u.
   !
   laue_expand_right = 100.0
/
ATOMIC_SPECIES
 Al  -1.0  Al.pbe-nl-rrkjus_psl.1.0.0.UPF
CELL_PARAMETERS angstrom
2.85671139599365 0.00000000000000 0.00000000000000
1.42835569799683 2.47398464021101 0.00000000000000
0.00000000000000 0.00000000000000 35.0000000000000
ATOMIC_POSITIONS angstrom
Al 1.4283556980 0.8246615467 -2.3324950875
Al 0.0000000000 1.6493230935 -0.0000000000
Al 0.0000000000 0.0000000000  2.3324950875
K_POINTS {automatic}
 8 8 1 1 1 0
SOLVENTS {mol/L}
 H2O  -1.0  H2O.spc.MOL
 Na+   1.0  Na+.aq.MOL
 Cl-   1.0  Cl-.aq.MOL

ESM-RISMの実行後、静電ポテンシャルと溶媒の電荷密度の結果を描画します。
溶質と溶媒の静電ポテンシャルは、左のパネルに示すように、真空領域と溶媒領域で傾きをもちます。しかし、全ポテンシャルは真空領域と溶媒領域で平坦になります。この結果は、表面電荷が、電解質のイオンに誘起された反対符号の電荷によって完全に遮蔽されていることを示しています。右のパネルでは、表面近くのCl\(^-\) の電荷密度が、電荷中性の表面の場合と比べて増加しています。一方、Na\(^+\) の電荷密度の変化は、Cl\(^-\) と逆の傾向を示します。これらの電荷密度分布の変化は、誘起された表面電荷と溶媒イオンの間のクーロン相互作用によるものです。

NaCl水溶液中の正に帯電したAl(111)スラブの静電ポテンシャルと電荷密度

2-3. NaCl水溶液中の負に帯電したAl(111)スラブに対するESM-RISM計算¶

3つ目として、\(\sigma=+0.02e\) の場合と比較するために、\(\sigma=-0.02e\) のAl(111)表面とNaCl水溶液の界面を計算します。

&control
   calculation  = 'scf'
   prefix       = 'Al111+0.02e.NaCl_aq'
   pseudo_dir   = './pseudo'
   outdir       = './output'
   trism        = .true.
/
&system
   ibrav = 0
   ntyp = 1, nat = 3
   ecutwfc = 25.0, ecutrho = 200.0
   nbnd = 16
   occupations = 'smearing', smearing = 'gauss', degauss = 0.01
   assume_isolated = "esm", esm_bc = "bc1", tot_charge = -0.02
/
&electrons
   electron_maxstep = 200
   mixing_beta = 0.1
   diagonalization = 'rmm', diago_rmm_conv = .false., diago_rmm_ndim = 4
/
&rism
   nsolv    = 3, closure  = 'kh'
   tempv    = 300.0  ! Kelvin
   ecutsolv = 160.0  ! Rydberg
   !
   ! Lennard-Jones for each atom
   !
   ! For Al using universal force field
   solute_lj(1) = 'uff'
   !
   ! 1D-RISM's setting
   !
   starting1d      = 'zero'
   rism1d_conv_thr = 1.0e-8
   rism1d_maxstep  = 10000
   mdiis1d_size    = 20
   mdiis1d_step    = 0.5
   !
   ! 3D-RISM's setting
   !
   starting3d      = 'zero'
   rism3d_maxstep  = 10000
   rism3d_conv_thr = 1.0e-6
   !
   ! For Laue-RISM calculation
   ! in a.u.
   !
   laue_expand_right = 100.0
/
ATOMIC_SPECIES
 Al  -1.0  Al.pbe-nl-rrkjus_psl.1.0.0.UPF
CELL_PARAMETERS angstrom
2.85671139599365 0.00000000000000 0.00000000000000
1.42835569799683 2.47398464021101 0.00000000000000
0.00000000000000 0.00000000000000 35.0000000000000
ATOMIC_POSITIONS angstrom
Al 1.4283556980 0.8246615467 -2.3324950875
Al 0.0000000000 1.6493230935 -0.0000000000
Al 0.0000000000 0.0000000000  2.3324950875
K_POINTS {automatic}
 8 8 1 1 1 0
SOLVENTS {mol/L}
 H2O  -1.0  H2O.spc.MOL
 Na+   1.0  Na+.aq.MOL
 Cl-   1.0  Cl-.aq.MOL

帯電した表面と電荷中性の表面のポテンシャルと密度の比較¶

電荷ゼロ

Al(111)スラブのポテンシャルと密度：電荷ゼロ

+0.02e

Al(111)スラブのポテンシャルと密度：電荷 +0.02e

−0.02e

Al(111)スラブのポテンシャルと密度：電荷 −0.02e

2-4. NaCl水溶液中で溶媒イオン（Na⁺）を陽に扱ったAl(111)スラブに対するESM-RISM計算¶

&control
   calculation  = 'scf'
   prefix       = 'Al111.NaCl_aq'
   pseudo_dir   = './pseudo'
   outdir       = './output'
   trism        = .true.
/
&system
!   ibrav = 1, a = 10.0
   ibrav = 0
   ntyp = 2, nat = 28
   ecutwfc = 40.0, ecutrho = 320.0
   occupations = 'smearing', smearing = 'gauss', degauss = 0.01
   assume_isolated = "esm", esm_bc = "bc1", tot_charge = 1.0
/
&electrons
   electron_maxstep = 200
   mixing_beta = 0.1
   diagonalization = 'rmm', diago_rmm_conv = .false., diago_rmm_ndim = 4
/
&rism
   nsolv    = 3
   closure  = 'kh'
   tempv    = 300.0  ! Kelvin
   ecutsolv = 160.0  ! Rydberg
   !
   ! Lennard-Jones for each atom
   !
   ! For Al using universal force field
   solute_lj(1) = 'uff'
   ! For Aqueous Na+
   solute_lj(2) = 'none'
   solute_epsilon(2) = 0.1301
   solute_sigma(  2) = 2.3500
   !
   ! 1D-RISM's setting
   !
   starting1d      = 'zero'
   rism1d_conv_thr = 1.0e-8
   rism1d_maxstep  = 10000
   mdiis1d_size    = 20
   mdiis1d_step    = 0.5
   !
   ! 3D-RISM's setting
   !
   starting3d      = 'zero'
   rism3d_maxstep  = 10000
   rism3d_conv_thr = 1.0e-6
   rism3d_conv_level = 1.0
   !
   ! For Laue-RISM calculation
   ! in a.u.
   !
   laue_expand_right = 100.0
/
ATOMIC_SPECIES
 Al  -1.0  Al.pbe-nl-rrkjus_psl.1.0.0.UPF 
 Na  -1.0  Na.pbe-spnl-rrkjus_psl.1.0.0.UPF

CELL_PARAMETERS angstrom
8.57013418798096 0.00000000000000 0.00000000000000
4.28506709399048 7.42195392063303 0.00000000000000
0.00000000000000 0.00000000000000 35.0000000000000

ATOMIC_POSITIONS angstrom
Al 1.4283556980 0.8246615467 -2.3324950875 
Al 4.2850670940 0.8246615467 -2.3324950875 
Al 7.1417784900 0.8246615467 -2.3324950875 
Al 2.8567113960 3.2986461869 -2.3324950875 
Al 5.7134227920 3.2986461869 -2.3324950875 
Al 8.5701341880 3.2986461869 -2.3324950875 
Al 4.2850670940 5.7726308272 -2.3324950875 
Al 7.1417784900 5.7726308272 -2.3324950875 
Al 9.9984898860 5.7726308272 -2.3324950875 
Al 0.0000000000 1.6493230935 0.0000000000 
Al 2.8567113960 1.6493230935 0.0000000000 
Al 5.7134227920 1.6493230935 0.0000000000 
Al 1.4283556980 4.1233077337 0.0000000000 
Al 4.2850670940 4.1233077337 0.0000000000 
Al 7.1417784900 4.1233077337 0.0000000000 
Al 2.8567113960 6.5972923739 0.0000000000 
Al 5.7134227920 6.5972923739 0.0000000000 
Al 8.5701341880 6.5972923739 0.0000000000 
Al 0.0000000000 0.0000000000 2.3324950875 
Al 2.8567113960 0.0000000000 2.3324950875 
Al 5.7134227920 0.0000000000 2.3324950875 
Al 1.4283556980 2.4739846402 2.3324950875 
Al 4.2850670940 2.4739846402 2.3324950875 
Al 7.1417784900 2.4739846402 2.3324950875 
Al 2.8567113960 4.9479692804 2.3324950875 
Al 5.7134227920 4.9479692804 2.3324950875 
Al 8.5701341880 4.9479692804 2.3324950875 
Na 6.4276006410 3.7109769603 12.3324950875 

K_POINTS {automatic}
2 2 1 1 1 0

SOLVENTS {mol/L}
 H2O  -1.0  H2O.spc.MOL
 Na+   1.0  Na+.aq.MOL
 Cl-   1.0  Cl-.aq.MOL

NaCl水溶液中で溶媒イオン（Na+）を陽に扱ったAl(111)スラブのポテンシャルと密度

2-5. NaCl水溶液中で溶媒イオン（Cl⁻）を陽に扱ったAl(111)スラブに対するESM-RISM計算¶

&control
   calculation  = 'scf'
   prefix       = 'Al111.NaCl_aq'
   pseudo_dir   = './pseudo'
   outdir       = './output'
   trism        = .true.
/
&system
   ibrav = 0
   ntyp = 2, nat = 28
   ecutwfc = 40.0, ecutrho = 320.0
   occupations = 'smearing', smearing = 'gauss', degauss = 0.01
   assume_isolated = "esm", esm_bc = "bc1", tot_charge = -1.0
/
&electrons
   electron_maxstep = 200
   mixing_beta = 0.1
   diagonalization = 'rmm', diago_rmm_conv = .false., diago_rmm_ndim = 4
/
&rism
   nsolv    = 3, closure  = 'kh'
   tempv    = 300.0  ! Kelvin
   ecutsolv = 160.0  ! Rydberg
   !
   ! Lennard-Jones for each atom
   !
   ! For Al using universal force field
   solute_lj(1) = 'uff'
   ! For Aqueous Cl-
   solute_lj(2) = 'none'
   solute_epsilon(2) = 0.1001
   solute_sigma(  2) = 4.4000
   !
   ! 1D-RISM's setting
   !
   starting1d      = 'zero'
   rism1d_conv_thr = 1.0e-8
   rism1d_maxstep  = 20000
   mdiis1d_size    = 20
   mdiis1d_step    = 0.5
   !
   ! 3D-RISM's setting
   !
   starting3d      = 'zero'
   rism3d_maxstep  = 10000
   rism3d_conv_thr = 1.0e-6
   rism3d_conv_level = 1.0
   !
   ! For Laue-RISM calculation
   ! in a.u.
   !
   laue_expand_right = 100.0
/
ATOMIC_SPECIES
 Al  -1.0  Al.pbe-nl-rrkjus_psl.1.0.0.UPF
 Cl  -1.0  Cl.pbe-nl-rrkjus_psl.1.0.0.UPF

CELL_PARAMETERS angstrom
8.57013418798096 0.00000000000000 0.00000000000000
4.28506709399048 7.42195392063303 0.00000000000000
0.00000000000000 0.00000000000000 35.0000000000000

ATOMIC_POSITIONS angstrom
Al 1.4283556980 0.8246615467 -2.3324950875
Al 4.2850670940 0.8246615467 -2.3324950875
Al 7.1417784900 0.8246615467 -2.3324950875
Al 2.8567113960 3.2986461869 -2.3324950875
Al 5.7134227920 3.2986461869 -2.3324950875
Al 8.5701341880 3.2986461869 -2.3324950875
Al 4.2850670940 5.7726308272 -2.3324950875
Al 7.1417784900 5.7726308272 -2.3324950875
Al 9.9984898860 5.7726308272 -2.3324950875
Al 0.0000000000 1.6493230935 0.0000000000
Al 2.8567113960 1.6493230935 0.0000000000
Al 5.7134227920 1.6493230935 0.0000000000
Al 1.4283556980 4.1233077337 0.0000000000
Al 4.2850670940 4.1233077337 0.0000000000
Al 7.1417784900 4.1233077337 0.0000000000
Al 2.8567113960 6.5972923739 0.0000000000
Al 5.7134227920 6.5972923739 0.0000000000
Al 8.5701341880 6.5972923739 0.0000000000
Al 0.0000000000 0.0000000000 2.3324950875
Al 2.8567113960 0.0000000000 2.3324950875
Al 5.7134227920 0.0000000000 2.3324950875
Al 1.4283556980 2.4739846402 2.3324950875
Al 4.2850670940 2.4739846402 2.3324950875
Al 7.1417784900 2.4739846402 2.3324950875
Al 2.8567113960 4.9479692804 2.3324950875
Al 5.7134227920 4.9479692804 2.3324950875
Al 8.5701341880 4.9479692804 2.3324950875
Cl 6.4276006410 3.7109769603 12.3324950875

K_POINTS {automatic}
2 2 1 1 1 0

SOLVENTS {mol/L}
 H2O  -1.0  H2O.spc.MOL
 Na+   1.0  Na+.aq.MOL
 Cl-   1.0  Cl-.aq.MOL

NaCl水溶液中で溶媒イオン（Cl-）を陽に扱ったAl(111)スラブのポテンシャルと密度

3. ESM-RISM計算による標準水素電極（SHE）電位の決定¶

グランドポテンシャルの計算¶

次の化学反応の左辺と右辺について、グランドポテンシャル（\(\Omega\)）のプロファイルを描きます。

\[ 1/2 \mathrm{H}_2 (\mathrm{gas}) +\mathrm{H}_2\mathrm{O} (\mathrm{S}) \rightleftharpoons \mathrm{H}_3\mathrm{O}^+ + e^- (\mathrm{M}) \]

左辺と右辺の状態の \(\Omega\)（\(\Omega_L\) と \(\Omega_R\)）の定義は次のとおりです。

\[ \begin{aligned} \Omega_{\mathrm{L}}&=1/2[A(\mathrm{H}_2, \mathrm{gas})+E_{\rm zp}(\mathrm{H}_2)]+\Omega(\mathrm{M}+\mathrm{H}_2\mathrm{O}, \mathrm{S})+E_{\rm zp}(\mathrm{H}_2\mathrm{O}) \\ \Omega_{\mathrm{R}}&=\Omega(\mathrm{M}+\mathrm{H}_3\mathrm{O}^+, \mathrm{S})+E_{\rm zp}(\mathrm{H}_3\mathrm{O}^+) \\ E(\mathrm{H}^+/\mathrm{H}_2)[\mathrm{vs. \ } E^{\rm S}_{\rm vac}]&=-\mu^{\mathrm{SHE}}_{\mathrm{e}} /\mathrm{e} [\mathrm{vs} \Phi_{\rm S}]. \end{aligned} \]

ここで、M は金属電極を表します。このチュートリアルでは、電極としてPt(111)表面を用いました。

溶媒分子を陽に扱った場合と扱わない場合のグランドポテンシャルのプロファイル¶

ここでは、溶媒（H\(_2\)OまたはH\(_3\)O\(^+\)）を陽に扱って電極電位を計算します。
入力ファイルを以下に示します。

入力ファイル：Pt(111)の清浄表面¶

tot_charge を-1.0から+0.8まで変化させました。

&control  
   calculation  = 'scf'  
   prefix       = '0.00'  
   pseudo_dir   = './pseudo'
   outdir       = './data'  
   trism        = .true.  
/  
&system  
   ibrav = 0  
   ntyp = 1, nat = 64   
   ecutwfc = 40.0, ecutrho = 320.0  
   occupations = 'smearing', smearing = 'gauss', degauss = 0.01  
   assume_isolated = "esm", esm_bc = "bc1", tot_charge = 0.0  
/  
&electrons  
   electron_maxstep = 200    
   mixing_beta = 0.1  
   diagonalization = 'rmm', diago_rmm_conv = .false., diago_rmm_ndim = 4  
/  
&rism  
   nsolv    = 3  
   closure  = 'kh'  
   tempv    = 300.0  ! Kelvin   
   ecutsolv = 160.0  ! Rydberg  
   !  
   ! Lennard-Jones for each atom  
   !  
   ! For Pt using fitted LJ-potential  
   !  
   solute_lj(1) = 'none',solute_epsilon(1) = 1.66, solute_sigma(1)=2.65   
   !  
   ! 1D-RISM's setting  
   !  
   starting1d      = 'zero'  
   rism1d_conv_thr = 1.0e-8  
   rism1d_maxstep  = 100000  
   mdiis1d_size    = 20  
   mdiis1d_step    = 0.1  
   !  
   ! 3D-RISM's setting  
   !  
   starting3d      = 'zero'  
   rism3d_maxstep  = 10000  
   rism3d_conv_thr = 1.0e-6  
   rism3d_conv_level = 0.5  
   mdiis3d_size    = 20  
   mdiis3d_step    = 0.8  
   ! 
   ! For Laue-RISM calculation
   ! in a.u.
   !
   laue_expand_right = 60.0
   laue_starting_right = -8.20
   laue_buffer_right = 10.0
/
ATOMIC_SPECIES
 Pt  -1.0  Pt.pbe-n-van.UPF

CELL_PARAMETERS angstrom
11.08743432900506 0.00000000000000 0.00000000000000
5.54371716450253 9.60199979171006 0.00000000000000
0.00000000000000 0.00000000000000 40.000

ATOMIC_POSITIONS angstrom
Pt 0.0000000000 0.0000000000 -10.0000000000 
Pt 2.7718585823 0.0000000000 -10.0000000000 
Pt 5.5437171645 0.0000000000 -10.0000000000 
Pt 8.3155757468 0.0000000000 -10.0000000000 
Pt 1.3859292911 2.4004999479 -10.0000000000 
Pt 4.1577878734 2.4004999479 -10.0000000000 
Pt 6.9296464556 2.4004999479 -10.0000000000 
Pt 9.7015050379 2.4004999479 -10.0000000000 
Pt 2.7718585823 4.8009998959 -10.0000000000 
Pt 5.5437171645 4.8009998959 -10.0000000000 
Pt 8.3155757468 4.8009998959 -10.0000000000 
Pt 11.0874343290 4.8009998959 -10.00000000000
Pt 4.1577878734 7.2014998438 -10.0000000000 
Pt 6.9296464556 7.2014998438 -10.0000000000 
Pt 9.7015050379 7.2014998438 -10.0000000000 
Pt 12.4733636201 7.2014998438 -10.00000000000
Pt 1.3859292911 0.8001666493 -7.7367869448  
Pt 4.1577878734 0.8001666493 -7.7367869448  
Pt 6.9296464556 0.8001666493 -7.7367869448  
Pt 9.7015050379 0.8001666493 -7.7367869448  
Pt 2.7718585823 3.2006665972 -7.7367869448  
Pt 5.5437171645 3.2006665972 -7.7367869448  
Pt 8.3155757468 3.2006665972 -7.7367869448  
Pt 11.0874343290 3.2006665972 -7.7367869448 
Pt 4.1577878734 5.6011665452 -7.7367869448  
Pt 6.9296464556 5.6011665452 -7.7367869448  
Pt 9.7015050379 5.6011665452 -7.7367869448  
Pt 12.4733636201 5.6011665452 -7.7367869448 
Pt 5.5437171645 8.0016664931 -7.7367869448  
Pt 8.3155757468 8.0016664931 -7.7367869448  
Pt 11.0874343290 8.0016664931 -7.7367869448 
Pt 13.8592929113 8.0016664931 -7.7367869448 
Pt -0.0000000000 1.6003332986 -5.4735738896 
Pt 2.7718585823 1.6003332986 -5.4735738896  
Pt 5.5437171645 1.6003332986 -5.4735738896  
Pt 8.3155757468 1.6003332986 -5.4735738896  
Pt 1.3859292911 4.0008332465 -5.4735738896  
Pt 4.1577878734 4.0008332465 -5.4735738896  
Pt 6.9296464556 4.0008332465 -5.4735738896  
Pt 9.7015050379 4.0008332465 -5.4735738896  
Pt 2.7718585823 6.4013331945 -5.4735738896  
Pt 5.5437171645 6.4013331945 -5.4735738896  
Pt 8.3155757468 6.4013331945 -5.4735738896  
Pt 11.0874343290 6.4013331945 -5.4735738896 
Pt 4.1577878734 8.8018331424 -5.4735738896  
Pt 6.9296464556 8.8018331424 -5.4735738896  
Pt 9.7015050379 8.8018331424 -5.4735738896  
Pt 12.4733636201 8.8018331424 -5.4735738896 
Pt 0.0000000000 0.0000000000 -3.2103608344  
Pt 2.7718585823 0.0000000000 -3.2103608344  
Pt 5.5437171645 0.0000000000 -3.2103608344  
Pt 8.3155757468 0.0000000000 -3.2103608344  
Pt 1.3859292911 2.4004999479 -3.2103608344  
Pt 4.1577878734 2.4004999479 -3.2103608344  
Pt 6.9296464556 2.4004999479 -3.2103608344  
Pt 9.7015050379 2.4004999479 -3.2103608344  
Pt 2.7718585823 4.8009998959 -3.2103608344  
Pt 5.5437171645 4.8009998959 -3.2103608344  
Pt 8.3155757468 4.8009998959 -3.2103608344  
Pt 11.0874343290 4.8009998959 -3.2103608344 
Pt 4.1577878734 7.2014998438 -3.2103608344  
Pt 6.9296464556 7.2014998438 -3.2103608344  
Pt 9.7015050379 7.2014998438 -3.2103608344  
Pt 12.4733636201 7.2014998438 -3.2103608344 

K_POINTS {automatic}
2 2 1 1 1 0

SOLVENTS {mol/L}
H2O  -1.0  H2O.tip5p.MOL
H3O+  1.0  H3O+.chuev.MOL
Cl-   1.0  Cl-.aq.MOL

入力ファイル：H₂O分子を置いたPt(111)¶

tot_charge を-1.0から0.0まで変化させました。

&control
   calculation  = 'scf'
   prefix       = '0.00'
   pseudo_dir   = './pseudo'
   outdir       = './data'
   trism        = .true.
/
&system
   ibrav = 0
   ntyp = 3, nat = 67
   ecutwfc = 40.0, ecutrho = 320.0
   occupations = 'smearing', smearing = 'gauss', degauss = 0.01
   assume_isolated = "esm", esm_bc = "bc1", tot_charge = 0.0
/
&electrons
   electron_maxstep = 200
   mixing_beta = 0.1
   diagonalization = 'rmm', diago_rmm_conv = .false., diago_rmm_ndim = 4
/
&rism
   nsolv    = 3
   closure  = 'kh'
   tempv    = 300.0  ! Kelvin
   ecutsolv = 160.0  ! Rydberg
   !
   ! Lennard-Jones for each atom
   !
   ! For Pt using fitted LJ-potential
   !
   solute_lj(1) = 'none',solute_epsilon(1) = 1.66, solute_sigma(1)=2.65 
   !
   ! For H
   !
   solute_lj(2) = 'none'
   solute_epsilon(2) = 0.046
   solute_sigma(  2) = 1.000
   !
   ! For O
   !
   solute_lj(3) = 'none'
   solute_epsilon(3) = 0.1600
   solute_sigma(  3) = 3.1200
   !
   ! 1D-RISM's setting
   !
   starting1d      = 'zero'
   rism1d_conv_thr = 1.0e-8
   rism1d_maxstep  = 100000
   mdiis1d_size    = 20
   mdiis1d_step    = 0.1
   !
   ! 3D-RISM's setting
   !
   starting3d      = 'zero'
   rism3d_maxstep  = 10000
   rism3d_conv_thr = 1.0e-6
   rism3d_conv_level = 0.5
   mdiis3d_size    = 20
   mdiis3d_step    = 0.8
   !
   ! For Laue-RISM calculation
   ! in a.u.
   !
   laue_expand_right = 60.0
   laue_starting_right = -8.20
   laue_buffer_right = 10.0
/
ATOMIC_SPECIES
 Pt  -1.0  Pt.pbe-n-van.UPF
 H   -1.0  H.pbe-rrkjus_psl.0.1.UPF
 O   -1.0  O.pbe-rrkjus.UPF

CELL_PARAMETERS angstrom
11.08743432900506 0.00000000000000 0.00000000000000
5.54371716450253 9.60199979171006 0.00000000000000
0.00000000000000 0.00000000000000 40.000

ATOMIC_POSITIONS angstrom
Pt 0.0000000000 0.0000000000 -10.0000000000  0 0 0
Pt 2.7718585823 0.0000000000 -10.0000000000  0 0 0
Pt 5.5437171645 0.0000000000 -10.0000000000  0 0 0
Pt 8.3155757468 0.0000000000 -10.0000000000  0 0 0
Pt 1.3859292911 2.4004999479 -10.0000000000  0 0 0
Pt 4.1577878734 2.4004999479 -10.0000000000  0 0 0
Pt 6.9296464556 2.4004999479 -10.0000000000  0 0 0
Pt 9.7015050379 2.4004999479 -10.0000000000  0 0 0
Pt 2.7718585823 4.8009998959 -10.0000000000  0 0 0
Pt 5.5437171645 4.8009998959 -10.0000000000  0 0 0
Pt 8.3155757468 4.8009998959 -10.0000000000  0 0 0
Pt 11.0874343290 4.8009998959 -10.0000000000  0 0 0
Pt 4.1577878734 7.2014998438 -10.0000000000  0 0 0
Pt 6.9296464556 7.2014998438 -10.0000000000  0 0 0
Pt 9.7015050379 7.2014998438 -10.0000000000  0 0 0
Pt 12.4733636201 7.2014998438 -10.0000000000  0 0 0
Pt 1.3859292911 0.8001666493 -7.7367869448  0 0 0
Pt 4.1577878734 0.8001666493 -7.7367869448  0 0 0
Pt 6.9296464556 0.8001666493 -7.7367869448  0 0 0
Pt 9.7015050379 0.8001666493 -7.7367869448  0 0 0
Pt 2.7718585823 3.2006665972 -7.7367869448  0 0 0
Pt 5.5437171645 3.2006665972 -7.7367869448  0 0 0
Pt 8.3155757468 3.2006665972 -7.7367869448  0 0 0
Pt 11.0874343290 3.2006665972 -7.7367869448  0 0 0
Pt 4.1577878734 5.6011665452 -7.7367869448  0 0 0
Pt 6.9296464556 5.6011665452 -7.7367869448  0 0 0
Pt 9.7015050379 5.6011665452 -7.7367869448  0 0 0
Pt 12.4733636201 5.6011665452 -7.7367869448  0 0 0
Pt 5.5437171645 8.0016664931 -7.7367869448  0 0 0
Pt 8.3155757468 8.0016664931 -7.7367869448  0 0 0
Pt 11.0874343290 8.0016664931 -7.7367869448  0 0 0
Pt 13.8592929113 8.0016664931 -7.7367869448  0 0 0
Pt -0.0000000000 1.6003332986 -5.4735738896  0 0 0
Pt 2.7718585823 1.6003332986 -5.4735738896  0 0 0
Pt 5.5437171645 1.6003332986 -5.4735738896  0 0 0
Pt 8.3155757468 1.6003332986 -5.4735738896  0 0 0
Pt 1.3859292911 4.0008332465 -5.4735738896  0 0 0
Pt 4.1577878734 4.0008332465 -5.4735738896  0 0 0
Pt 6.9296464556 4.0008332465 -5.4735738896  0 0 0
Pt 9.7015050379 4.0008332465 -5.4735738896  0 0 0
Pt 2.7718585823 6.4013331945 -5.4735738896  0 0 0
Pt 5.5437171645 6.4013331945 -5.4735738896  0 0 0
Pt 8.3155757468 6.4013331945 -5.4735738896  0 0 0
Pt 11.0874343290 6.4013331945 -5.4735738896  0 0 0
Pt 4.1577878734 8.8018331424 -5.4735738896  0 0 0
Pt 6.9296464556 8.8018331424 -5.4735738896  0 0 0
Pt 9.7015050379 8.8018331424 -5.4735738896  0 0 0
Pt 12.4733636201 8.8018331424 -5.4735738896  0 0 0
Pt 0.0000000000 0.0000000000 -3.2103608344  0 0 0
Pt 2.7718585823 0.0000000000 -3.2103608344  0 0 0
Pt 5.5437171645 0.0000000000 -3.2103608344  0 0 0
Pt 8.3155757468 0.0000000000 -3.2103608344  0 0 0
Pt 1.3859292911 2.4004999479 -3.2103608344  0 0 0
Pt 4.1577878734 2.4004999479 -3.2103608344  0 0 0
Pt 6.9296464556 2.4004999479 -3.2103608344  0 0 0
Pt 9.7015050379 2.4004999479 -3.2103608344  0 0 0
Pt 2.7718585823 4.8009998959 -3.2103608344  0 0 0
Pt 5.5437171645 4.8009998959 -3.2103608344  0 0 0
Pt 8.3155757468 4.8009998959 -3.2103608344  0 0 0
Pt 11.0874343290 4.8009998959 -3.2103608344  0 0 0
Pt 4.1577878734 7.2014998438 -3.2103608344  0 0 0
Pt 6.9296464556 7.2014998438 -3.2103608344  0 0 0
Pt 9.7015050379 7.2014998438 -3.2103608344  0 0 0
Pt 12.4733636201 7.2014998438 -3.2103608344  0 0 0
H 8.3155757468 4.8009998959 6.7896391656 
H 6.8046757468 4.8009998959 6.7896391656 
O 7.5601257468 5.3899498959 6.7896391656

K_POINTS {automatic}
2 2 1 1 1 0

SOLVENTS {mol/L}
H2O  -1.0  H2O.tip5p.MOL
H3O+  1.0  H3O+.chuev.MOL
Cl-   1.0  Cl-.aq.MOL

入力ファイル：H₃O⁺イオンを置いたPt(111)¶

tot_charge を0.0から1.0まで変化させました。

&control
   calculation  = 'scf'
   prefix       = '0.00'
   pseudo_dir   = './pseudo'
   outdir       = './data'
   trism        = .true.
/
&system
   ibrav = 0
   ntyp = 3, nat = 68
   ecutwfc = 40.0, ecutrho = 320.0
   occupations = 'smearing', smearing = 'gauss', degauss = 0.01
   assume_isolated = "esm", esm_bc = "bc1", tot_charge = 1.0
   starting_charge(2)=0.33
/
&electrons
   electron_maxstep = 200
   mixing_beta = 0.1
   diagonalization = 'rmm', diago_rmm_conv = .false., diago_rmm_ndim = 4
/
&rism
   nsolv    = 3
   closure  = 'kh'
   tempv    = 300.0  ! Kelvin
   ecutsolv = 160.0  ! Rydberg
   !
   ! Lennard-Jones for each atom
   !
   ! For Pt using fitted LJ-potential
   !
   solute_lj(1) = 'none',solute_epsilon(1) = 1.66, solute_sigma(1)=2.65
   !
   ! For H
   !
   solute_lj(2) = 'none'
   solute_epsilon(2) = 0.046
   solute_sigma(  2) = 1.000
   !
   ! For O
   !
   solute_lj(3) = 'none'
   solute_epsilon(3) = 0.1600
   solute_sigma(  3) = 3.1200
   !
   ! 1D-RISM's setting
   !
   starting1d      = 'zero'
   rism1d_conv_thr = 1.0e-8
   rism1d_maxstep  = 100000
   mdiis1d_size    = 20
   mdiis1d_step    = 0.1
   !
   ! 3D-RISM's setting
   !
   starting3d      = 'zero'
   rism3d_maxstep  = 10000
   rism3d_conv_thr = 1.0e-6
   rism3d_conv_level = 0.5
   mdiis3d_size    = 20
   mdiis3d_step    = 0.8
   !
   ! For Laue-RISM calculation
   ! in a.u.
   !
   laue_expand_right = 60.0
   laue_starting_right = -8.20
   laue_buffer_right = 10.0
/
ATOMIC_SPECIES
 Pt  -1.0  Pt.pbe-n-van.UPF
 H   -1.0  H.pbe-rrkjus_psl.0.1.UPF
 O   -1.0  O.pbe-rrkjus.UPF

CELL_PARAMETERS angstrom
11.08743432900506 0.00000000000000 0.00000000000000
5.54371716450253 9.60199979171006 0.00000000000000
0.00000000000000 0.00000000000000 40.000

ATOMIC_POSITIONS angstrom
Pt 0.0000000000 0.0000000000 -10.0000000000  0 0 0
Pt 2.7718585823 0.0000000000 -10.0000000000  0 0 0
Pt 5.5437171645 0.0000000000 -10.0000000000  0 0 0
Pt 8.3155757468 0.0000000000 -10.0000000000  0 0 0
Pt 1.3859292911 2.4004999479 -10.0000000000  0 0 0
Pt 4.1577878734 2.4004999479 -10.0000000000  0 0 0
Pt 6.9296464556 2.4004999479 -10.0000000000  0 0 0
Pt 9.7015050379 2.4004999479 -10.0000000000  0 0 0
Pt 2.7718585823 4.8009998959 -10.0000000000  0 0 0
Pt 5.5437171645 4.8009998959 -10.0000000000  0 0 0
Pt 8.3155757468 4.8009998959 -10.0000000000  0 0 0
Pt 11.0874343290 4.8009998959 -10.0000000000  0 0 0
Pt 4.1577878734 7.2014998438 -10.0000000000  0 0 0
Pt 6.9296464556 7.2014998438 -10.0000000000  0 0 0
Pt 9.7015050379 7.2014998438 -10.0000000000  0 0 0
Pt 12.4733636201 7.2014998438 -10.0000000000  0 0 0
Pt 1.3859292911 0.8001666493 -7.7367869448  0 0 0
Pt 4.1577878734 0.8001666493 -7.7367869448  0 0 0
Pt 6.9296464556 0.8001666493 -7.7367869448  0 0 0
Pt 9.7015050379 0.8001666493 -7.7367869448  0 0 0
Pt 2.7718585823 3.2006665972 -7.7367869448  0 0 0
Pt 5.5437171645 3.2006665972 -7.7367869448  0 0 0
Pt 8.3155757468 3.2006665972 -7.7367869448  0 0 0
Pt 11.0874343290 3.2006665972 -7.7367869448  0 0 0
Pt 4.1577878734 5.6011665452 -7.7367869448  0 0 0
Pt 6.9296464556 5.6011665452 -7.7367869448  0 0 0
Pt 9.7015050379 5.6011665452 -7.7367869448  0 0 0
Pt 12.4733636201 5.6011665452 -7.7367869448  0 0 0
Pt 5.5437171645 8.0016664931 -7.7367869448  0 0 0
Pt 8.3155757468 8.0016664931 -7.7367869448  0 0 0
Pt 11.0874343290 8.0016664931 -7.7367869448  0 0 0
Pt 13.8592929113 8.0016664931 -7.7367869448  0 0 0
Pt -0.0000000000 1.6003332986 -5.4735738896  0 0 0
Pt 2.7718585823 1.6003332986 -5.4735738896  0 0 0
Pt 5.5437171645 1.6003332986 -5.4735738896  0 0 0
Pt 8.3155757468 1.6003332986 -5.4735738896  0 0 0
Pt 1.3859292911 4.0008332465 -5.4735738896  0 0 0
Pt 4.1577878734 4.0008332465 -5.4735738896  0 0 0
Pt 6.9296464556 4.0008332465 -5.4735738896  0 0 0
Pt 9.7015050379 4.0008332465 -5.4735738896  0 0 0
Pt 2.7718585823 6.4013331945 -5.4735738896  0 0 0
Pt 5.5437171645 6.4013331945 -5.4735738896  0 0 0
Pt 8.3155757468 6.4013331945 -5.4735738896  0 0 0
Pt 11.0874343290 6.4013331945 -5.4735738896  0 0 0
Pt 4.1577878734 8.8018331424 -5.4735738896  0 0 0
Pt 6.9296464556 8.8018331424 -5.4735738896  0 0 0
Pt 9.7015050379 8.8018331424 -5.4735738896  0 0 0
Pt 12.4733636201 8.8018331424 -5.4735738896  0 0 0
Pt 0.0000000000 0.0000000000 -3.2103608344  0 0 0
Pt 2.7718585823 0.0000000000 -3.2103608344  0 0 0
Pt 5.5437171645 0.0000000000 -3.2103608344  0 0 0
Pt 8.3155757468 0.0000000000 -3.2103608344  0 0 0
Pt 1.3859292911 2.4004999479 -3.2103608344  0 0 0
Pt 4.1577878734 2.4004999479 -3.2103608344  0 0 0
Pt 6.9296464556 2.4004999479 -3.2103608344  0 0 0
Pt 9.7015050379 2.4004999479 -3.2103608344  0 0 0
Pt 2.7718585823 4.8009998959 -3.2103608344  0 0 0
Pt 5.5437171645 4.8009998959 -3.2103608344  0 0 0
Pt 8.3155757468 4.8009998959 -3.2103608344  0 0 0
Pt 11.0874343290 4.8009998959 -3.2103608344  0 0 0
Pt 4.1577878734 7.2014998438 -3.2103608344  0 0 0
Pt 6.9296464556 7.2014998438 -3.2103608344  0 0 0
Pt 9.7015050379 7.2014998438 -3.2103608344  0 0 0
Pt 12.4733636201 7.2014998438 -3.2103608344  0 0 0
H 8.6355757468 5.7429998959 6.7896391656
H 8.6355757468 4.3299998959 7.6056391656
H 8.6355757468 4.3299998959 5.9736391656
O 8.3155757468 4.8009998959 6.7896391656  0 0 0

K_POINTS {automatic}
2 2 1 1 1 0

SOLVENTS {mol/L}
H2O  -1.0  H2O.tip5p.MOL
H3O+  1.0  H3O+.chuev.MOL
Cl-   1.0  Cl-.aq.MOL

グランドポテンシャルのプロファイルの結果¶

溶媒を陽に扱った場合（左）と扱わない場合（右）のグランドポテンシャルのプロファイル

左：溶媒を陽に扱った場合の \(\Omega(\mu)\)
右：溶媒を陽に扱わない場合の \(\Omega(\mu)\)

左右のパネルの交点は、それぞれ \(\mu=-5.23\)eVと \(\mu=-5.26\)eVです。

電極を用いないSHE電位¶

SHE電位は、Helmholtz自由エネルギーの差からも計算できます。

\[ \begin{aligned} G_{\mathrm{L}}&=1/2[A(\mathrm{H}_2, \mathrm{gas})+E_{\rm zp}(\mathrm{H}_2)]+A(\mathrm{H}_2\mathrm{O}, \mathrm{S})+E_{\rm zp}(\mathrm{H}_2\mathrm{O}) \\ G_{\mathrm{R}}&=A(\mathrm{H}_3\mathrm{O}^+, \mathrm{S})+E_{\rm zp}(\mathrm{H}_3\mathrm{O}^+) \end{aligned} \]

結果は次のとおりです。

Ry2eV=27.2116*0.5
TS_H2=-0.404
A_H2=-2.3322957476*Ry2eV+TS_H2
A_H2O=-34.2914608381*Ry2eV
A_H3O=-35.1019053889*Ry2eV

Ezp_H2=0.270
Ezp_H2O=0.562
Ezp_H3O=0.909

G_L=0.5*(A_H2+Ezp_H2)+A_H2O+Ezp_H2O
G_R=A_H3O+Ezp_H3O

print('E(H+/H2)[vs. E_vac] (eV) = '+str(-(G_L-G_R)))

E(H+/H2)[vs. E_vac] (eV) = 5.2536282720733425

この値（E(H+/H2)[vs. E_vac]=5.25eV）は、グランドポテンシャルのプロファイルから得た交点とほぼ一致します。
入力ファイルも以下に示します。

入力ファイル：H₂O分子¶

&CONTROL
  calculation  = 'relax'       ! 'scf' 'relax'
  outdir       = './data'
  pseudo_dir   = './pseudo'
  prefix       = 'H2O.HCl_aq'
  trism        = .true.
!  lfcp         = .true.
  max_seconds  =  72000
  nstep        = 100
/
&SYSTEM
  ibrav = 1, A = 20.0
  ntyp = 2, nat = 3
  nbnd = 20
  ecutwfc = 40, ecutrho = 320
  occupations = "smearing", smearing = "gauss", degauss = 0.01
  assume_isolated = "esm", esm_bc = "bc1", tot_charge = 0.0
!  starting_charge(2) = 0.33
/
&ELECTRONS
  electron_maxstep = 500, mixing_beta = 0.1, !conv_thr = 1.0d-8
  diagonalization = 'rmm', diago_rmm_conv = .false., diago_rmm_ndim = 4
!  startingwfc = 'file'
!  startingpot = 'file'
/
&IONS
  ion_dynamics = 'bfgs'
!  upscale = 1.d0
/
&RISM
  nsolv = 3, closure  = 'kh'
  tempv = 300, ecutsolv = 160
!Lennard-Jones
  solute_lj(1) = 'none', solute_epsilon(1) = 0.1554, solute_sigma(1) = 3.166 !TIP5P O
  solute_lj(2) = 'none', solute_epsilon(2) = 0.0460, solute_sigma(2) = 1.000 !modified H
!Laue-RISM
  laue_expand_left = 60
  laue_expand_right = 60
  rism3d_maxstep = 10000
  rism3d_conv_thr = 1.0d-6
  rism3d_conv_level = 0.5
!  starting1d = "file"
!  starting3d = "file"
/

K_POINTS gamma

SOLVENTS {mol/L}
H2O  -1.0  H2O.tip5p.MOL
H3O+  1.0  H3O+.chuev.MOL
Cl-   1.0  Cl-.aq.MOL

ATOMIC_SPECIES
O   16.00  O.pbe-rrkjus.UPF 
H   1.008  H.pbe-rrkjus_psl.0.1.UPF

ATOMIC_POSITIONS {angstrom}
O        0.   0.   0.
H        0.6  0.6  0.6
H        0.6  0.6 -0.6

入力ファイル：H₃O⁺イオン¶

&CONTROL
  calculation  = 'relax'       ! 'scf' 'relax'
  outdir       = './data'
  pseudo_dir   = './pseudo'
  prefix       = 'H2O.HCl_aq'
  trism        = .true.
!  lfcp         = .true.
  max_seconds  =  72000
  nstep        = 100
/
&SYSTEM
  ibrav = 1, A = 20.0
  ntyp = 2, nat = 4
  nbnd = 20
  ecutwfc = 40, ecutrho = 320
  occupations = "smearing", smearing = "gauss", degauss = 0.01
  assume_isolated = "esm", esm_bc = "bc1", tot_charge = 1.0
!  starting_charge(2) = 0.33
/
&ELECTRONS
  electron_maxstep = 500, mixing_beta = 0.1, !conv_thr = 1.0d-8
  diagonalization = 'rmm', diago_rmm_conv = .false., diago_rmm_ndim = 4
!  startingwfc = 'file'
!  startingpot = 'file'
/
&IONS
  ion_dynamics = 'bfgs'
!  upscale = 1.d0
/
&RISM
  nsolv = 3, closure  = 'kh'
  tempv = 300, ecutsolv = 160
!Lennard-Jones
  solute_lj(1) = 'none', solute_epsilon(1) = 0.1554, solute_sigma(1) = 3.166 !TIP5P O
  solute_lj(2) = 'none', solute_epsilon(2) = 0.0460, solute_sigma(2) = 1.000 !modified H
!Laue-RISM
  laue_expand_left = 60
  laue_expand_right = 60
  rism3d_maxstep = 10000
  rism3d_conv_thr = 1.0d-6
  rism3d_conv_level = 0.5
!  starting1d = "file"
!  starting3d = "file"
/

K_POINTS gamma

SOLVENTS {mol/L}
H2O  -1.0  H2O.tip5p.MOL
H3O+  1.0  H3O+.chuev.MOL
Cl-   1.0  Cl-.aq.MOL

ATOMIC_SPECIES
O   16.00  O.pbe-rrkjus.UPF 
H   1.008  H.pbe-rrkjus_psl.0.1.UPF

ATOMIC_POSITIONS {angstrom}
O        0.   0.   0.
H        0.6  0.6  0.6
H        0.6  0.6 -0.6
H       -0.6 -0.5  0.