溶液中の電子移動反応の古典的分子動力学（MD) 計算:定式化 Formulation ofClassical Molecular Dynamics Simulation of Electron Transfer Reaction in Solution

山本雅博
Modiﬁed on 平成 18 年 6 月 7 日 5 : 39 pm（以下は山本がこれまで報告されている電子移動の分子動力学計算の論文を読んで Self-study したものをまとめたものです。未保証ですので，間違っていたらごめんなさい。）
水溶液中で，酸化還元を行うイオンは水に溶媒和されており，その自由エネルギーは溶媒分子の構造的な（距離・配向）揺らぎにともない，時間とともに変化している。このイオンが何らかの原因（電極からの電子移動，バルク中での酸化還元反応等）により電子を授受してその価数が変わるとする。
この電子移動は，溶媒の揺らぎ（ps オーダー）に比べて早い反応（おそらく fs レベル）であり，従って，電子移動反応の反応経路は，溶媒の揺らぎを止めた状態でのイオンの電荷の変化に対応すると考えてよい。古典的な分子動力学計算では，イオンと溶媒の相互作用は，１）静電相互作用と２）van derWaals 相互作用（Lennard-Jones ポテンシャル）の和になる。
酸化体と還元体で２）の相互作用は同じにとるので，１）の相互作用の大小により速い電子移動の反応経路は決定される。
今，１電子酸化還元反応を Oxz + e ⇀（溶↽ Redz−1 とし，媒和環境を固定した状況下における）反応座標に相当する ∆Eを以下のように定義する。

∆E= −(Vox,sol − VRed,sol )∑ (z − 1)zje2 ∑ zzj
−= −4πϵ0|rj − R|
|rj − R|jje ∑ zj e
= −4πϵ0 j |rj − R|!
O( kB T 2 CV )図 1: ∆E の時間変化の模式図。∆E = ∆e になる確率分布より自由エネルギー F (∆e) が求まる。

し，バイアスのかかったハミルトニアンを H ′ = H + Vb とする。ある物理量の熱平均は，∫dxAe−βH∫〈A〉H =(5)dxe−βH∫βV ∫dxAe−β(H+Vb )e b dxe−β(H+Vb )∫∫=(6)dxe−βHdxe−β(H+Vb )∫βV ∫dxAe−β(H+Vb )e b dxe−βVb e−βH∫∫=(7)dxe−βHdxe−β(H+Vb )(1)= 〈AeβVb 〉H ′ 〈e−βVb 〉H∆E は時間ともに揺らぐ (rj (t))。十分に長い分子動力学計
算を行えば，∆E がある値 ∆e をとる確率で,P (∆e) を求めることができる。(図１参照）
P (∆e) = 〈δ(∆e − ∆E)〉で与えられる。
A = δ(∆e − ∆E) とし Vb = Vb (∆E) とすると，(2)〈δ(∆e − ∆E)〉H

ここで，δ(x) は Dirac のデルタ関数である。系の自由エネル
ギーを F (∆e) とすると，確率 P (∆e) は，

P (∆e) = exp[−βF (∆e)],

β=1kB T= 〈δ(∆e − ∆E)eβVb (∆E) 〉H ′ 〈e−βVb 〉H= eβVb (∆e) 〈δ(∆e − ∆E)〉H ′ 〈e−βVb 〉H
(9)

(3)

ここで，デルタ関数の性質を使い，H ′ における熱平均には無関係な項 eβVb (∆e) を，くくり出した。自由エネルギーは，以下の様に与えられる。

(4)F (∆e) = −kB T ln〈δ(∆e − ∆E)〉H= −kB T ln〈δ(∆e − ∆E)〉H ′ − Vb (∆e) + C(10)で与えられる。従って，F (∆e) = −kB T ln〈δ(∆e − ∆E)〉(8)となる。長時間の分子動力学計算を行い，上の式から自由エネルギーを求まることができる。ただし，∆E の揺らぎの幅は，√kB T のオーダーである。(N V T アンサンブルの場合 kB T 2 CVのオーダーである。水の場合 CV =75 J K−1 mol−1 at 300 Kとすると 300 K でのエネルギーの揺らぎは，77 meV = 7.5
kJ mol−1 となる。）
多くの場合，電子移動の活性化障壁は，kB T のオーダーよりも大きく，通常の MD をおこなっていても，めったに訪れないエネルギー領域にある。このために，バイアスとなるポテンシャルを系に強制的にかけて，特別の領域の ∆e に系を置き，その効果を後で差し引くことにより F (∆e) を求めることが できる。(umbrella sampling) 今，系のバイアスのかかってない Hamiltonian（古典的には全エネルギー）を H とC は，Vb の選択によって決まる定数である。Vb としては，２次関数が良く用いられる。Vb (∆E) = (k/2)(∆E −∆e0 )2 . k, ∆e0を調整して，多数の”overlapping windows(or bin)”で計算を行いそれぞれをつなぎあわせて，バイアスのかかってないなめら
かな F (∆e) を求める。（WHAM(Weighted Histogram Analysis Method) という方法で決定できる。）
バイアスがかかった時の，MD 計算では通常の分子間力に加えて以下の力も考慮する必要がある。i 番目の原子における力を Fi とすると，
H′1= H+k(∆E − ∆e0 )22(11)Fi∇ri∑jzj e|rj − R|= −∇ri H ′(12)= −∇ri H − k(∆E − ∆e0 )∇ri ∆E= −∇ri H∑ zj e−e−k(∆E − ∆e0 )∇ri 4πϵ0 j |R − rj |= (i∂∂∂+j+k)∂xi∂yi∂zi∑(13)(zj e)/j[(xj − Rx )2 + (yj − Ry )2 + (zj − Rz )2 ]1/2
= zi e[(xi − Rx )i + (yi − Ry )j + (zi − Rz )k][(xi − Rx )2 + (yi − Ry )2 + (zi − Rz )2 ]−3/2
Vb = (k/2)(rI − r0I )2 または Vb = (k/2)(zI − zI0 )2 とす
ればイオンの位置に”overlapping windows(or bin)”をかけた
umbrella sampling により自由エネルギーの位置依存性が計算できる。
H′FIH′FI

図 4: Upper panel: solid lines give the MD-generated freeenergy surfaces for Cl2− , Cl− , Cl, and Cl+ , dashed lines
give the parabolic ﬁts of the minima. Lower panel: average number of hydrogen bonds of a water molecule in theﬁrst solvation shell as a function of the generalized solvent
coordinate. J. Chem. Phys. 115, pp. 8540–8546, 2001 C.
Hartnig and M. T. M. Koperk(rI − r0I )22= −∇rI H ′ = −∇rI H − k(rI − r0I )k= H + (zI − zI0 )22= −∇rI H ′ = −∇rI H − k(zI − zI0 )k= H+(14)

図 2: Macus 理論の模式図。横軸の反応座標は，MD 計算の場合エネルギーの単位になる。

図 5: (a) The Cl–O pair radial distribution function as afunction of the generalized solvent coordinate. Numbers indicate the corresponding values of e (in kJ/mol). (b) The
location of the ﬁrst maximum in the Cl–O pair radial distribution function as a function of the generalized solvent
coordinate. Numbers in the lower part, which correspondto the charge of the ion (atom), indicate the location of theminimum in the free-energy curves. J. Chem. Phys. 115,
pp. 8540–8546, 2001 C. Hartnig and M. T. M. Koper.

図 3: 反応座標 x(= ∆e) での確率分布およびそれから計算した自由エネルギーの模式図。J. Chem. Phys. 1990 93, pp.
8682–8692.G. King and A. Warshel2