Classical Molecular Dynamics Simulation of Metal Electrodes-Electrolyte Interface

Abstract

電極電解液界面ではバルクの溶液とは大きく異なる特異な溶媒構造のもとで多彩かつ重要な反応が進行する．古典動力学シミュレーションは反応場としての電極電解液界面に対する微視的な理解を得ることが出来る有効なアプローチの一つである．本稿では，金属電極電解質水溶液系の古典動力学シミュレーションを行うにあたり注意すべき三点: 相互作用ポテンシャルの選択，電極の分極と定電位条件の記述，長距離静電相互作用の評価 について述べる．

1 はじめに

現代の科学技術や工学に欠かせない存在となっているリチウムイオン電池や 次世代のナトリウムイオン電池は，電極電解液界面で起こる反応によって 性能が大きく左右される．この反応を理解するためには，界面における酸化還元体と，電極や周囲の溶媒分子，イオンとの静電相互作用を 考慮することが重要である．溶媒と電極が特異な相互作用をすることで，溶媒分子が界面に緩く束縛された特有の溶媒構造が現れる．酸化還元体は，バルク溶液とは大きく異なるこの非等方的な溶媒構造の中で，複雑な溶媒効果を受けながら電極と電子をやり取りする．さらに，電圧を印加すると静電場により界面溶媒構造が変化するため，酸化還元体が受ける溶媒効果は複雑に変化すると予想されるが，これらの詳細は十分には明らかになっていない．このような基礎的な問題をミクロなレベルで解明することは，学術的意義は勿論のこと，新世代の優れた電池を開発する上でも重要である．

電子のやり取りに対する界面の影響を見積もるには，界面における静電ポテンシャルの変化を見るのが良い．古くから用いられるアプローチの一つで，溶媒を連続誘電体として扱うGouy-Chapman-Sternモデルを適用すると，静電ポテンシャルはFigure 1のように電極表面からバルクに向かって単調に変化する． [1]

Figure 1.

Electrostatic potenial profile obtained from Gouy-Chapman-Stern model for an electrode-electrolyte interface.

一方，溶媒分子を顕わに扱って評価した静電ポテンシャルは，界面特有の溶媒構造を反映し，Figure 2のように複雑に変化する．電極と酸化還元体間の電子のやり取りは，この静電ポテンシャルが大きく変化する辺りで起こりやすいと考えられるため，そのキネティクスを正しく理解するためには溶媒分子を顕わに考慮した界面の溶媒構造を知る必要がある．

Figure 2.

Electrostatic potenial profile obtained from an all atom simulation.

電極電解液界面に対する電子・原子レベルの理解を得ることは，AFMや分光といった界面に適用出来る強力な実験手法を用いても 難しく，原子レベルのレゾリューションを持つ理論計算からのアプローチが重要な役割を果たす．例えば分子動力学法のようなシミュレーションによって，原子レベルでの界面溶媒構造，溶媒効果，その電極表面からの距離依存性，またそれらに対する印加電圧の影響を調べることが出来る．このとき満たすべき条件として，

1. 溶媒の熱揺らぎを考慮すること(配置サンプリングを十分に行うこと)

2. 原子間，特に電極-電解液間の相互作用を正確に記述すること

3. 非等方的な系の記述，特に長距離クーロン相互作用の適切な扱いが挙げられる．最初の条件についてここで少し詳しく述べておきたい．統計サンプリングの問題は溶媒が関係する全シミュレーションにおいて重要であり，バルク溶液系の研究では注意深く議論されることが多い．電極電解液界面ではバルク溶液に比べ溶媒構造の緩和が遅く，十分なサンプリングはより困難である．

Figure 3はPt電極-水界面において，Na⁺イオンを電極に近づけるときの 自由エネルギー変化(平均力ポテンシャル)を古典分子動力学計算により評価したものである．Na⁺イオン-電極間距離が約2.5 Åにおける平均力ポテンシャルの値が，用いるトラジェクトリの長さ(10 ∼500 ns)によって10 kJ/mol 以上変わってしまうことが見て取れる．また十分に収束した結果を得るためには，400 ns以上の長さのトラジェクトリが必要であることもわかる．

Figure 3.

Potential of mean force profiles for a sodium cation approaching a negatively charged Pt electrode in aqueous solution calculated from trajectories ranging from 10 ns to 500 ns. The horizontal axis shows the distance from the negatively charged electrode.

本稿では，上記の条件を満たしうる方法である 電極電解液界面に対する古典動力学法に関して，我々が行った研究 [2,3] から得た経験を参考にしながら，バルク溶液系にはない問題とそれに対するアプローチについて述べる．系として意識するのは金属-水界面に絞るが，電極や電解液が他のものに変わった場合も共通点は少なくないだろう．

まず次節では，界面での化学過程を考えるにあたり基礎的かつ重要な金属表面近傍の水の構造について，実験と第一原理計算からわかっていることを簡単に述べる．次に古典分子動力学計算を行うにあたり，特徴的な界面溶媒構造を再現するような 金属電極-水分子間の相互作用ポテンシャルについて述べる．その後，金属電極を扱う上で重要と考えられている金属の分極効果と，エネルギー準位だけでなく溶媒構造も変化させる印加電圧の効果を，電圧一定条件の下でどのように取り入れるかに対するアプローチについて議論する．更に電極電解液系に対してよく用いられるスラブモデルを扱うにあたり，長距離クーロン相互作用をEwald法で評価する際の注意点について述べる．

2 金属表面近傍の水の構造

この節では主にHodgsonとHaqらのレビュー [4] とSakongらの第一原理計算結果 [5,6]に基づき，金属表面上の水の構造について簡単にまとめる．常温での電極電解液系の界面溶媒構造を実験的に調べるのは難しいが，低温における金属表面への水数分子の吸着と単分子層は，走査型トンネル顕微鏡 (STM)，X線光電子分光法 (XPS)，振動スペクトル，低速電子線回折 (LEED)，ヘリウム原子散乱といった実験的手法と，密度汎関数理論 (DFT) に基づく計算によって調べられている．まずどの金属が用いられるか，また表面にどの格子面が用いられるかで，反応性や界面溶媒構造が大きく異なりうることを意識する必要がある．例えば実験や計算でよく調べられているPt(111)表面は安定で反応性が低く，印加電圧が低く吸着酸素がないならば水分子はそのまま吸着する [7]. 反応性の低い金属，またその低指数面では，吸着エネルギーと水分子間の水素結合エネルギーが 大きく違わないため，それらの繊細なバランスによって界面構造が決まることになる． 例えばPt (111)面だと水一分子の吸着エネルギーが0.35 eV程度 [8]で， 対して気相中で水素結合した水二分子の結合エネルギーが0.22 eVである [9]. また水素結合した水分子間の酸素原子間距離が表面金属原子間距離と大体同じであることも，特徴的な溶媒構造(二次元の水素結合ネットワーク)が界面で出来上がることに関係している．水一分子だと金属原子の真上に酸素原子が来るように吸着する.水分子は分子面が金属表面と大体平行になり，水素原子を金属表面とは反対側に少し向けた構造をとる．これは水分子の孤立電子対を金属原子の空軌道と 共有するためと理解されている．水六分子が(111)面に吸着すると，全ての酸素原子が金属原子上に位置したヘキサマーが形成される．水単分子層になると，ヘキサマーのときと同様に酸素原子が金属原子上に位置し二次元に拡がった構造を形成するか，或いはより良い層内水素結合ネットワークの形成が優先され必ずしも酸素原子が金属原子上に位置しない構造が出来る，ともいわれている．またこのとき水素原子が金属表面側を向くか外側を向くかについても議論がある．常温の電解液の場合は，より多くの水分子との相互作用と熱揺らぎも 界面構造に影響を与えるため，状況は一層複雑になる．以上のことから，シミュレーションによって意味のある議論をする為には，水分子間に加えて水分子と金属表面間の相互作用を適切に記述しなければならない．

有限温度における電極電解液界面の溶媒構造を調べるための有効な手法として 第一原理分子動力学計算(AIMD)が用いられる．原子数が多い系に対し溶媒配置の サンプリングを行うため，計算コストが比較的小さくて済むDFTが用いられることが殆どである．このとき汎関数に何を用いるかや分散力を別に考慮するかが結果に大きく影響するため，それらの慎重な選択が重要になる．汎関数として例えばRPBE [10]を用い，GrimmeのD3分散力補正 [11] を加えると，CCSD (T)レベルでの水分子クラスター安定構造とエネルギー相対値，バルクの水の動径分布関数，拡散係数の実験値を良く再現することが示されている．金属表面への水一分子の吸着エネルギーも計算され，汎関数依存性，分散力補正の有無の影響が大きいことがわかっているが，対応する実験値が無くどの計算レベルが適当かは明らかになっていない． [6] AIMD計算により界面溶媒分布関数が評価されてはいるが， [5,7] その莫大な計算コストのためにトラジェクトリの長さは数10 ps程度に限られており，界面の溶媒構造について決定的な結果を得るのは容易ではない．

3 金属-溶媒間相互作用ポテンシャル

古典分子動力学法は電極電解液界面を調べるために必要な大量の統計サンプリングを可能とし，常温での界面溶媒構造やその溶媒効果を調べることの出来るアプローチである．但し信頼出来る第一原理計算により評価された電解液間，電極-電解液間の相互作用を定量的に再現出来ることが必要である．ここでは特にPt電極表面-水分子間相互作用の経験ポテンシャルについて述べたい．第2節で述べたように，水一分子の安定吸着構造では，金属原子の真上(on-topサイト)に酸素原子が位置し，水分子面がPt表面とほぼ平行になる．これを経験力場で記述するとき，通常の力場のようにクーロン以外の相互作用の記述に LJポテンシャルを用いてしまうと，on-topサイトではなく金属原子間の上(hollowサイト，或いはbridgeサイト) に酸素原子が位置してしまう.これは二体のLJポテンシャルでは近接原子の数が多くなるような位置で安定化されるためである [12] ．このため，金属原子-水分子間におけるクーロン以外の特異な相互作用を記述するための 専用の相互作用ポテンシャルを構築する必要がある．

このようなポテンシャルの例として，まずSiepmannとSprikによるものが挙げられる [12]. これは電極表面原子と酸素原子間の距離だけでなく，水分子の双極子モーメントが金属原子と酸素原子を結ぶベクトルとのなす角度も用い，さらに金属二原子-酸素原子の三体ポテンシャルを導入することで，金属クラスターと水分子に対する量子化学計算結果を 再現するよう作成されたものである．このポテンシャルを用いてPt(111)表面上に水一分子を吸着させると，酸素原子がon-topサイトに位置するが水素原子が表面から離れすぎる(水分子の双極子モーメントが表面に垂直な方向を向く)ように なっており，MichaelidesらのPW91レベルの計算 [8]や SakongらのRPBE-D3レベルの計算結果 [6]とは一致しない．この原因は，フィットに用いた量子化学計算レベル(extended Hückel)の違いが考えられる．高精度な量子化学計算結果を再現するようにパラメタをフィットし直すことで，より信頼性の高い結果が得られると期待できる．

SiepmannとSprikのポテンシャルは最近でもよく使われている [13, 14]. ここではこのポテンシャルを用いた古典分子動力学計算結果の一例として，Pt電極と1M K₃Fe (CN)₆電解質水溶液からなる電極電解液系について 得られた溶媒分布関数をFigure 4に示す．z = 0, 50 Åの位置に表面がある二つのPt電極間には2 Vの電圧をかけている．表面から約2.7 Åの位置に水の吸着層に由来する酸素原子のピークが鋭く現れている．負に帯電した右側の電極から約5 Å離れた位置にK${}^{+}$イオンのピークがあり，この位置にK${}^{+}$が多く存在することがわかる．また表面から約7 Å離れた位置にFe (CN) ₆^3-の中心原子Feが多く存在する．この分布関数は異なる初期配置を持つ1 nsのトラジェクトリ10本から平均を取って得たものである．

Figure 4.

Number density profiles along the direction normal to the surfaces of the Pt metal electrodes. The system is composed of a positively charged Pt electrode (left electrode), a negatively charged Pt electrode (right electrode), and 1 M K₃Fe(CN)₆ aqueous solution, with 2 V voltage applied between the two electrodes. The profiles for atoms Fe in Fe(CN)₆^3-, K, O in water molecules are shown.

その他のポテンシャル関数として，水分子の双極子モーメントと金属表面のなす角度を取り入れたもの， [15] on-top, hollow, bridgeサイトの違いを表現するため金属表面原子の並びを表す余弦関数を取り入れたもの [16] などが挙げられる．また特異な金属原子-溶媒分子間相互作用の記述に加え，溶媒分子の分解反応も表せるようにした ポテンシャル， [17, 18] 大量の構造に対する第一原理エネルギー計算結果を 正確に再現でき，溶媒分子の分解も表せるニューラルネットワークポテンシャルなども作成されている [19].

4 金属電極の分極と定電位条件

極性をもった溶媒分子やイオンが作りだす電場に応答し，電極が電子分極して表面に電荷が誘起される．誘起された電荷は界面の溶媒構造や界面で起こる種々の反応に影響を与えるため重要であり， [20,21] 電極電解液界面のモデル化において取り入れられることが多い．また電極電解液系の実際の状況と合わせるため，電極間電位差が一定に保たれるよう工夫が必要である．電位差は電極電荷と電解液分子・イオンの分布によって決まることから，電解液の原子配置が時々刻々と変化するのに応答し，電位差一定の条件の下で電極電荷も変化することになる．アプローチとしては，古典電磁気学の鏡像法に基づき鏡像電荷を導入するものと，金属内部の静電ポテンシャルが 一定であるとして金属誘起原子電荷を決定するものがある.後に述べるように，後者は溶液系の分子動力学計算で 分子の電子分極を表現するためにしばしば用いられる化学ポテンシャル平衡法と考え方は同じである．金属電極の記述の仕方として，電極原子を顕わに考えずに完全導体平板とみなすものと，電極原子を顕わに扱い，分極を表現するため金属原子電荷を導入するものがある．電極原子の並びが界面溶媒構造に大きな影響を与えるので，界面の詳細な解析を行いたい場合は 金属原子を顕わに考慮したアプローチが望ましい.更に後者は平行平板電極だけでなく 任意の形状の電極を表し易いことも利点として挙げられる．

前者では，平板表面に平行なx, y方向には周期的，表面に垂直なz方向には非周期的であることを境界条件を通して取り入れたGreen関数を用いて，金属電極が作りだす静電ポテンシャルを評価する アプローチが挙げられる [22,23]. 印加電圧も境界条件を調整することで取り入れられる．Petersenらは金属原子を顕わに扱いつつ，鏡像法に基づく少し簡単化された方法を提案している [24].

後者では，電極原子電荷を可変とし，電極原子上の静電ポテンシャルが各電極内で一定という条件から 電荷を決定するアプローチがよく用いられる.二つある電極にかかる静電ポテンシャルの差が印加電圧となる．SiepmannとSprikは，電極原子電荷を点電荷ではなく，拡がりξをパラメタに持つ電荷q_iのGauss型の分布関数   

${\rho }_i\left(r\right)=\frac{q_i}{{\left(2\pi {\xi }^2\right)}^{3/2}}exp\left[-\frac{{(r-r_i)}^2}{2{\xi }^2}\right]$(1)

で表し，金属電極内のエネルギー変化がGauss型原子電荷間のクーロン相互作用(遮蔽クーロン相互作用)と 自己相互作用とからなるとした [12]. このときの金属電極内の静電エネルギーは次式で表される．   

$U=\frac{1}{2}\sum _{i\ne j}{q}_i{q}_j\frac{\text{erf}\left(r_{ij}/2\xi \right)}{r_{ij}}+\frac{1}{2}\sum _i\frac{q_i^2}{\sqrt{\pi }\xi }$(2)

また金属原子にかかる静電ポテンシャルとポテンシャル一定条件は次のように表せる．   

$\sum _{j(\ne i)}{q}_j\frac{\text{erf}\left(r_{ij}/2\xi \right)}{r_{ij}}+\frac{q_i}{\sqrt{\pi }\xi }=V_{\text{const}}$(3)

左辺第二項は，点電荷ではなく拡がりを持つ電荷分布の自己相互作用に由来する．

この方法は分子系の分極可能力場で用いられる化学ポテンシャル平衡法 (chemical potential equalization method，或いはelectronegativity equalization method, またはfluctuating charge model) [25, 26] を金属電極に適用したもの [27, 28]と捉えることも出来る．化学ポテンシャル平衡法では，エネルギーを原子電荷の変化{q_i}で二次まで展開した次の式で表す．   

$U=U_0+\sum _i{\chi }_0{q}_i+\frac{1}{2}\sum _{i,j}{q}_i{J}_{ij}{q}_j$(4)

χ₀は原子電荷が変化する前の電気陰性度，J_ijは対角項(i = j)が"硬さ"(chemical hardness)，非対角項(i ≠ j)は遮蔽クーロン相互作用に対応する．原子電荷でエネルギーを微分して得られる各原子の化学ポテンシャルが 電極内で等しくなるとする:   

$({\chi }_0+\sum _{j(\ne i)}{J}_{ij}{q}_j)+J_{ii}{q}_i=-{\mu }_{\text{const}}.$(5)

原子電荷として式(1)を採用し，エネルギーを静電エネルギーのみに限ると，J_ijは 式(2)の遮蔽クーロン相互作用関数   

$J_{ij}=\frac{\text{erf}(r_{ij}/2\xi )}{r_{ij}},$(6)

J_iiは自己相互作用関数   

$J_{ii}=\frac{1}{\sqrt{\pi }\xi }$(7)

に一致する．式(5)がSiepmannとSprikの方法における静電ポテンシャル一定条件の式(3)に対応する．印加電圧の効果は二つの電極間の化学ポテンシャルの差が印加電圧に等しくなるよう，化学ポテンシャルの値を適切にシフトすればよい．或いは電気陰性度 χ₀をシフトすると考えることも出来る．すると式(4)右辺第二項がReedらの方法 [29, 30] における"electrical work term"に対応することがわかる．

変化する金属原子電荷の値は静電ポテンシャル，或いは化学ポテンシャル一定の連立方程式 (3)または(5)を解く，または拡張ラグランジアンを導入し電荷の運動方程式を評価すれば得ることが出来る [12]. 完全導体としてモデル化された金属電極では，電荷の大部分が表面第一層に誘起されるため，三層ほどでよくモデル化出来る．以上のアプローチを用いることで，金属電極の電子分極の効果を記述出来ると同時に，一定の電圧の下でのシミュレーションが可能となる．

Figure 5は，2 Vの印加電圧をかけたPt 電極と 1 M K₃Fe (CN)₆電解質水溶液からなる電極電解液系 (Figure 4と同じ系)における静電ポテンシャル変化 (ポアソンポテンシャルプロファイル)を評価したものである．二つの電極表面 (z = 0.5 Å) 間で 確かに電位差が印加電圧に等しくなっていることがわかる．電極表面から約10 Å離れた辺りまで特異な界面溶媒構造のため電位が振動し，溶媒を連続誘電体として考えるGouy-Chapman-Sternモデルから予想される単調な変化とは大きく異なる．このポアソンポテンシャルプロファイルは，まずSiepmannとSprikの方法を用いた定電位古典動力学計算を行い，Figure 4の溶媒分布関数に原子電荷をかけて得られる電荷分布関数ρ_q (z) を求め，金属電極表面電荷密度も含めた上で 一次元ポアソン方程式${\partial }^2/\partial z^2\varphi \left(z\right)=-{\rho }_q\left(z\right)/{ϵ}_0$ の積分形 (SI単位系)   

$\varphi \left(z\right)=\varphi \left(0\right)-\frac{1}{{ϵ}_0}{\int }_0^z{\rho }_q\left(z^{\prime }\right)(z-z^{\prime })d{z}^{\prime }$(8)

を数値的に評価することで得られる．なおFigure 5はFigure 4の計算と同じく，独立な10本の1 nsトラジェクトリを用いて得た．

Figure 5.

Poisson potential profile along the direction normal to the surfaces of the Pt metal electrodes calculated for the same system as described in Figure 4.

5 スラブモデルでの長距離クーロン相互作用の評価

極性溶媒とイオンからなる電解液，そして金属の分極を扱うことから，電極電解液界面を調べるためにはEwald法を用いて長距離静電相互作用を評価する必要がある．Figure 6に示すようなスラブモデルでは，電極表面に平行な方向(x, y方向)には無限に拡がっているとし，単純に周期境界条件を課して表現することが出来る．一方表面に垂直な方向(z方向)は非周期的であり同じ扱いは出来ない．このような系に対しては二次元Ewald法を用いるべきであるが， [29] 統計サンプリングをなるべく多量に行う為には，計算効率がより高い三次元Ewald法の使用が望まれる． [31,32,33] しかしその場合，z方向のイメージセルとの非物理的なクーロン相互作用が含まれてしまう．これは電極の外側に挿入した真空層を大きくとろうが取り除くことは難しい．この問題は印加電圧を課しz方向に大きな双極子モーメントが生じる場合に顕著に現れ，z方向に余計な電場がかかることで界面溶媒構造やポアソンポテンシャルに大きな影響を与えてしまう．

Figure 6.

 A slab model for studying an electrode-electrolyte interface. This model consists of two electrode layers (gold spheres), solvent (red and white spheres), and vacuum layers, all of which are contained in the primary cell (its boundaries are indicated with black lines).

この問題は，三次元Ewald法で表面項を考慮すれば回避出来る．簡単のため点電荷から構成される系を考える.表面項も含めた三次元Ewald法による静電エネルギーは   

$\begin{array}{c}U=\frac{1}{2}\sum _i\sum _j\sum _n{q^{\prime }}_i{q}_j\frac{\text{erfc}\left(\alpha \left|r_{ij}+nL\right|\right)}{\left|r_{ij}+nL\right|}\\ +\frac{1}{2\pi V}\sum _i\sum _j\sum _{k(\ne 0)}{q}_i{q}_j\frac{4{\pi }^2}{k^2}exp\left(-\frac{k^2}{4{\alpha }^2}\right)\\ \times cos(k\cdot r_{ij})\\ -\frac{\alpha }{\sqrt{\pi }}\sum _i{q}_i^2+S\left(M\right)\end{array}$(9)

で表さる．i, jは原子のラベル，n = (n_x, n_y, n_z)でn_x, n_y, n_zは整数である．右辺第一項目のプライムはn = 0のときi = jの項を除くことを意味する．また$nL=(n_{\text{x}}{L}_{\text{x}},n_{\text{y}}{L}_{\text{y}},n_{\text{z}}{L}_{\text{z}})$， $k=(2\pi m_{\text{x}}/L_{\text{x}},2\pi m_{\text{y}}/L_{\text{y}},2\pi m_{\text{z}}/L_{\text{z}})$で，Lx, Ly, Lzはそれぞれx, y, z方向のシミュレーションセルの長さ，m_x, m_y, m_zは整数である．右辺最後の項が表面項で基本セルが持つ双極子モーメントMに依存する． [31,34,35,37] バルク溶液系では表面項を省いてしまうことも多いが，非等方的な系で大きな双極子モーメントを持つような場合は考慮する必要がある．表面項の形は，Ewald法における本来は条件収束であるクーロン相互作用の無限級数和の処理の仕方で決まる．体積Vの基本セルを球状に複製するようにクーロン相互作用を足し合わせ，その外側が真空であるとするなら表面項は   

$S\left(M\right)=\frac{2\pi }{3V}{M}^2$(10)

となる．完全導体で囲まれているとするなら表面項はゼロとなり，これがいわゆるtin-foil boundary conditionである．基本セル複製の順序をスラブモデルの形状に対応させ，まずx, y平面上に複製し，その後z方向に複製していく場合，その外側を真空とするなら表面項は   

$S\left(M\right)=\frac{2\pi }{V}{M}_z^2$(11)

(M_zはMのz成分)となる． [35,38] また表面項は位置zにおける静電ポテンシャルに   

$V_{\text{surf}}\left(z\right)=-\frac{2\pi }{V}\sum _i{q}_i{(z-z_i)}^2$(12)

の寄与をする [36,37,38]. 電極の外側の真空層をある程度大きくとり (2枚の電極表面間距離程度 [31,39])，エネルギーに式(11)を加えることで，計算コストを殆ど増やすことなく，三次元Ewald法で二次元Ewald法とほぼ同じ 結果を得ることができ極めて有用である．YehとBerkowitz [39] が有効性を示したこの補正は，多くの電極電解液界面シミュレーションで用いられている．

6 まとめ

本稿では，バルク溶液にはない，電極電解液界面の古典動力学シミュレーションに 特有の三つの事柄について簡単にまとめた．電解液と接した電極の表面電子状態がどうなるのか，溶媒の化学吸着による電極-電解液間で電荷がどれ程移動するのか，溶媒の分解反応が起こりうるか，などを意識するなら，第一原理分子動力学計算を行うことが一つの理想といえる．しかし現時点では，それらと同等の重要性を持つといえる電解液の熱揺らぎの効果を 十分に取り入れることは極めて困難である．古典動力学シミュレーションでは電子が関係する情報の殆どが欠落してしまうが，適切な相互作用ポテンシャルの使用，電解液に対する電極の電子的応答(電子分極)の記述，長距離静電相互作用の評価を適切に行うことにより，熱揺らぎの効果を十分に含めた上で，複雑な界面の構造や反応場としての性質を調べることが出来る．十分な統計サンプリングが取れるという古典動力学法のメリットを活かすことで，界面に対する新たな知見がこれからも得られることだろう．

謝辞

本稿で紹介した計算結果は松三勇介氏の努力により得られたものである．本研究はJSPS科研費18K14179と，"京都大学 実験と理論計算科学のインタープレイによる触媒・電池の元素戦略研究拠点" (ESICB, Kyoto University)の助成を受けて行われた．

References

• [1]   Allen J Bard and Larry R Faulkner, Electrochemical methods: fundamentals and applications. New York: Wiley, 2 edition, 2001.
• [2]   Y. Matsumi, H. Nakano, H. Sato, Chem. Phys. Lett., 681, 80 (2017).
• [3]   H. Nakano, Y. Matsumi, H. Sato, in press.
• [4]   A. Hodgson, S. Haq, Surf. Sci. Rep., 64, 381 (2009).
• [5]   S. Sakong, K. Forster-Tonigold, A. Groß, J. Chem. Phys., 144, 194701 (2016).
• [6]   S. Sakong, M. Naderian, K. Mathew, R. G. Hennig, A. Groß, J. Chem. Phys., 142, 234107 (2015).
• [7]   J. Le, A. Cuesta, J. Cheng, J. Electroanal. Chem., 819, 87 (2018).
• [8]   A. Michaelides, V. A. Ranea, P. L. de Andres, D. A. King, Phys. Rev. Lett., 90, 216102 (2003).
• [9]   B. Temelso, K. A. Archer, G. C. Shields, J. Phys. Chem. A, 115, 12034 (2011).
• [10]   B. Hammer, L. B. Hansen, J. K. Nørskov, Phys. Rev. B, 59, 7413 (1999).
• [11]   S. Grimme, J. Antony, S. Ehrlich, H. Krieg, J. Chem. Phys., 132, 154104 (2010).
• [12]   J. I. Siepmann, M. Sprik, J. Chem. Phys., 102, 511 (1995).
• [13]   A. P. Willard, S. K. Reed, P. A. Madden, D. Chandler, Faraday Discuss., 141, 423 (2009).
• [14]   D. T. Limmer, A. P. Willard, P. A. Madden, D. Chandler, J. Phys. Chem. C, 119, 24016 (2015).
• [15]   S. N. Steinmann, R. Ferreira De Morais, A. W. Götz, P. Fleurat-Lessard, M. Iannuzzi, P. Sautet, C. Michel, J. Chem. Theory Comput., 14, 3238 (2018).
• [16]   K. Raghavan, K. Foster, K. Motakabbir, M. Berkowitz, J. Chem. Phys., 94, 2110 (1991).
• [17]   K.-Y. Yeh, M. J. Janik, J. K. Maranas, Electrochim. Acta, 101, 308 (2013).
• [18]   N. Onofrio, A. Strachan, J. Chem. Phys., 143, 054109 (2015).
• [19]   S. K. Natarajan, J. Behler, Phys. Chem. Chem. Phys., 18, 28704 (2016).
• [20]   A. A. Lee, S. Perkin, J. Phys. Chem. Lett., 7, 2753 (2016).
• [21]   J. Vatamanu, D. Bedrov, O. Borodin, Mol. Simul., 43, 838 (2017).
• [22]   K. Takae, A. Onuki, J. Phys. Chem. B, 119, 9377 (2015).
• [23]   A. P. dos Santos, M. Girotto, Y. Levin, J. Chem. Phys., 147, 184105 (2017).
• [24]   M. K. Petersen, R. Kumar, H. S. White, G. A. Voth, J. Phys. Chem. C, 116, 4903 (2012).
• [25]   A. K. Rappe, W. A. Goddard, J. Phys. Chem., 95, 3358 (1991).
• [26]   S. W. Rick, S. J. Stuart, B. J. Berne, J. Chem. Phys., 101, 6141 (1994).
• [27]   C. G. Guymon, R. L. Rowley, J. N. Harb, D. R. Wheeler, Condens. Matter Phys., 8, 335 (2005).
• [28]   T. Liang, A. C. Antony, S. A. Akhade, M. J. Janik, S. B. Sinnott, J. Phys. Chem. A, 122, 631 (2018).
• [29]   S. K. Reed, O. J. Lanning, P. A. Madden, J. Chem. Phys., 126, 084704 (2007).
• [30]   S. K. Reed, P. A. Madden, A. Papadopoulos, J. Chem. Phys., 128, 124701 (2008).
• [31]   A. Arnold, J. de Joannis, C. Holm, J. Chem. Phys., 117, 2496 (2002).
• [32]   S. Tyagi, A. Arnold, C. Holm, J. Chem. Phys., 129, 204102 (2008).
• [33]   A. Arnold, C. Holm, Comput. Phys. Commun., 148, 327 (2002).
• [34]   S. W. de Leeuw, J. W. Perram, E. R. Smith, Proc.- Royal Soc.A, Math. Phys. Eng. Sci., 373, 27 (1980).
• [35]   E. R. Smith, J. S. Rowlinson, Proc.- Royal Soc. A, Math. Phys. Eng. Sci., 375, 475 (1981). DOI:10.1098/rspa.1981.0064
• [36]   H. D. Herce, A. E. Garcia, T. Darden, J. Chem. Phys., 126, 124106 (2007).
• [37]   M. Mazars, Phys. Rep., 500, 43 (2011).
• [38]   A. P. dos Santos, M. Girotto, Y. Levin, J. Chem. Phys., 144, 144103 (2016).
• [39]   I. -C. Yeh, M. L. Berkowitz, J. Chem. Phys., 111, 3155 (1999).