1　はじめに
済液内や蛋白質中における酸化湿元中心間の電子移動反応を考える。電子移動を１次の非断熱遷移過程と仮定し、摂動論（Fermiの黄金則）とCondon近似（波動関数の電子成分と原子核成分の分離）を適用すると、反応速度定数は次式で表される。1)
H』は電子結合積分、FCはFranck-Condon因子にエネルギー準位密度を掛けて熱平均をとったものである。本稿では、この後者を分子計算によって求める方法について考察する。上式は、古典または高温極限において熱活性化型
に帰着する。△F°は、反応の自由エネルギー変化、入は再配置エネルギーである。これらは、反応の始伏態と終伏態のポテンシャル力場を決定すれば、モンテカルロ法や分子動力学法などにより計算することが出来る。さらに、後述の久保-豊沢による時間相関関数形式²」と分子動力学法を併用すれば、式(1)のFC因子そのものを量子論あるいは半古典論的に計算することも可能である。これらの分子計算において、媒質の電子分極効果を取り入れる有用な方法の一つが、本特集号のテーマである電荷平衡法である。我々の行った具体的なシミュレーション手法と結果については、第３節で紹介する。特に重要なのは、溶媒の電子分極効果によって再配置エネルギーが修正（renormalization）を受ける様子を、酸化湿元中心間距離(Rda)の関数として明らかにした点にある。これは、誘電連続体モデルにより導いかれた式1）
(3)における光学誘電率ε ∞の効果として定性的に理解されるが、微視的な分子モデルに基づく計算は従来存在しなかった。電荷平衡法がこの分野で有効に利用された一例である。また、adとaaは電子供与体と受容体の有効半径を表すが、これらの数値を例えば分子体積から見積もったとしても、誘竜体モデルでは定量性は期待できない。分光法によって入を見積もる試みもなされているが、測定が容易でない系がむしろ多数派であるようなので、やはり分子計算への期待は大きい。なお、ερは静的誘電率、Δeは移動電荷量である。 1これらの計算と結果について記述する前に、第２節において、電子移動反応における媒質の電子分極効果について若干の考察を加える。本稿で紹介する手法では、酸化湿元中心間の電子的結合と媒質電子分極とを独立に扱っているが、その妥当性に関しては注意を要するからである。ここでは、両者の相対的な時間スケールが関わっている。第３節で、具体的な計算手法を述べる。極性溶媒のように静電相互作用が強く揺らぎの大きい系の動力学計算においては、標準的な電荷平衡法は本来的な不安定性を争んでいる。ここでは、パラメータ数を増やさずにこれを回避する方法を紹介する。続いて、パラメータの決定手続き、分子シミュレーションによる再配置エネルギーの計算、久保-豊沢公式によるFC因子の計算について概説する。計算結果の一部を紹介した後、最終節にまとめと将来展望を簡単に記す。
2
電子移動反応における媒質の電子分極効果について
2.1 エネルギー準位配置と時間スケール
酸化湿元中心間距離の大きい長距離電子移動系では、媒質の空軌道を介した電子移動や占有軌道を介したホール移動が起こり得るので、電子移動過程と媒質の電子分極効果とは必ずしも独立に扱えない。しかし、次節以下で扱う水溶媒や、より一般に透明液体中では、媒質の電子励起エネルギーが可視光領域よりも大きいことから推察されるように、媒質を介した電子移動とホール移動が同時に大きく寄与する可能性は低い。特に、酸化湿元中心の電子準位が媒質のHOMO-LUMOギャップ内の中程にある場合には、電子移動過程と媒質電子分極を独立に扱う近似が妥当であると期待できる。そうでない場合には、本稿で記述するような電荷平衡法の単純な適用には注意を要する。以下では、上の仮定が成り立つとする。一般に、電子移動反応は分子構造の変化を伴うが、いま仮に反応の遷移状態で分子構造を凍結したとすると、酸化還元中心の電子状態は|H all/ftの時間周波数で振動する。実際には分子構造の緩和が起こり、それに伴い電子状態の局在が起こるので、この振動は減衰する(decoherence)が、 |H el|/hは電子移動の時間スケールの良い目安になる。これに対し、媒質の電子分極の時間スケールは、簡単に見積もるならば、その媒質の電子励起エネルギーをћω。とした場合のω。が対応するが、透明液体の場合にはこれは可視領域（数eV）を超える。一方、H elはおよそkcal/molのオーダーである。（H elがそれよりも大きい場合には、そもそも電子移動を1次の非断熱過程として捉える描像が破綻しかねない。）すなわち、ここで考えている様な系については、
frue ≫ | Hel |
(4)
が成り立つ。この場合、媒質の電子分極の時間スケールは電子移動過程よりも十分に速く、いわば後者に断熱的に追随すると見なすことが出来る。これに対し、逆の大小関係が成り立つ場合には、酸化湿元中心間の電子交換が速く、媒質の電子分極はその平均場に応答するような状況になる。後者は、電荷移動錯体を形成する場合などの強い相互作用系に相当する。次節以下では、上式が成り立つような弱結合系を取り扱うが、強結合系との違いを明らかにしておくことは、竜荷平衡法シミュレーションを適用する際にも重要であろう。この点については、第３節でも再度触れる。
2.2 ガウス型分極モデルによる経路積分解析
ここでは、前節で定性的に考察した状況について、簡単なモデルを設定し、なるべく厳密に取り扱ってみる。（読者の興味に依っては、飛ばして次節に進んで頂いても構いません。）まず、電子移動系を次の二状態ハミルトニアンで記述する。
反応の始状態と終状態を( i )と( j )で表し、m iとm fは対応する双極子結合係数、Δεは、媒質を除いた反応系に固有のエネルギー差、Jは前頁の|H allと同じものである。H aは媒質（熱浴）のハミルトニアン、Eは媒質からの電場で、これらは媒質の座標×（核と電子を含む）に依存する。この電場εを媒質の核座標成分、すなわち分子構造の揺らぎに起因する成分επと、電子分極成分ε。との和に分解する。さらに、熱浴ハミルトニアンについては、次の調和型モデルを仮定
化により、系の基子力学的分配関数を以下のように計算することができる。まず、式(5)は、次のスピン-ポゾン型に書き換えられる。
(7)ここで、H b ≡ H b + π E + Δε/2、Δm ≡ m i − m f、π ≡ (m i + m f )/2。σ xとσ zは通常のパウリ行列である。力学変数としてξπ(t)、ξρ(t)、および電子移動の状態を表すσ(t)（これは、±1の二値のみを取る）を考え、系の分配関数Q = Tr[e−βH]を虚数時間の経路積分表示で書くと、次のようになる。
(8)通常の経路積分法の通り、Sはβhの虚数時間でループ状に閉じた経路における作用である。 ƒ DE n (t)とƒ DE e (t)は連続変数の経路積分、∑ (g(t)}は±1の二値を断続的にスイッチするようなスピン経路の総和を表す。式(6)を代入し、ε o(t)の積分において２次のキュミュラント展開を取ると、上式は次のように変形できる。
Q。は、式(6)の第1項の電子分極部分εε/2α。のみをハミルトニアンとする仮想系の分配関数で、上式における統計平均(…)。も、これを参照したものである。5[σ(t),En(t)]は、電子分極Eo(t)に関する積分を先に実行したことにより修正(renormalization)を受けた作用である。これに含まれる核成分の分極率、双極子結合係数、エネルギーギャップも、α n ≡ α n /(1 − α eαnk 2 )、 Δm′ ≡ Δm(1 + kαe)、Δε′ ≡ Δε + αe π Δmのように修正を受ける。表式は省略するが、熱浴ハミルトニアンH Bも同様の修正を受け、それを以下F Bとする。つまり、電子分極E (t)に関する統計平均の効果は、♂を始めとする修正された諸量に繰り込まれ、電子分極の揺らぎを表す時間相関関数が、式(9)の最終指数因子に残った形になっている。ここまでの変形は厳密である。さて、式(4)に対応するような、速い電子分極の極限において、(Ee(t)Ee(t'))eをhaeδ(t − t')で近似することにする。これにより、ε (t)は式(9)から消え、分配関数は作用写によって決まることになる。すなわち、電子分極効果を繰り込んだ有効ハミルトニアンは、
となる。これは、式(7)と同型だが、全分極εが核成分εに置き換わり、Δm、Δεなども、ε eの効果を繰り込んだプライム( ′ )付きの量に換わっている。次節で説明する電荷平衡法シミュレーションでは、電荷分布は分子構造変化にほぼ断熱的に追随するとしている。すなわち、速い電子分極の極限である式(10)のモデルに対応していることになる。しかし、後述の運動方程式(15)の導出から分かるように、このシミュレーションにおける電荷自由度の運動は、正確な量子力学的時間発展とは無関係である。したがって、次節の電荷平衡法シミュレーションによって式(9)中のξ｡(t)の時間相関関数を議論することは出来ない。ともあれ、Δm′、Δε′などがα。を含んでいることに示されるように、式(10)の豆は電子分極効果を実効的に取り込んでいる。詳細は省くが、この互によって再配置エネルギーと自由エネルギー変化を計算すると、次のようになる。
α。＝ 0の極限を取ったものが、電子分極を無視(non-polarizable)したシミュレーションから得られる量に対応する。この場合の再配置エネルギー λnpは、式(11)でプライム(')を取り除いたものになる。ただし、電子分極を無視したシミュレーションでは、溶媒分子の双極子モーメントは気相よりも大きい値を用いる。これは、α nの値を補正することに相当する。この結果得られる再配置エネルギーは、電子分極を取り入れたときの値よりも大きくなる。一方、式(3)の誘電連続体モデルにおいては、電子分極を無視した再配置エネルギーはε∞ = 1
である。過去の文献では、同じ因子1/ϵ。で△F 0をもスケールするべきだという説も存在するがう)、我々の式(11)はこれを否定している。以上、式(5, 6)によるモデル的考察ではあったが、分子計算において電子分極効果を実効的に取り入れるということの意味について、定性的な理解を深める一助になれば幸いである。火節から、具体的なシミュレーションについて述べる。
電子移動反応の分子動力学シミュレーション
3
3.1 動的電荷平衡法における不安定性と修正モデル
本特集号で再三説明されているように、電荷平衡法では、系のポテンシャル関数を原子座標の組{r}および点電荷の組{Q}によって次式で表す。(簡単のため、点電荷は原子上にあるとする。)
第２項はクーロン相互作用、第３項はそれ以外の原子間相互作用ポテンシャルである。後者にはレナード・ジョーンズ関数がよく用いられる。第１項が電荷の関数としての原子のポテンシャルエネルギーであり、通常は次の２次式で表される。
χ⁰とη⁰は、原子に固有の電気陰性度と硬さ(hardness)を表す定数である。分子系における電気陰性度と硬さは、χ a = ∂U/∂Q a、2η a = ∂²U/∂Q aで与えられる。これらは分子構造や化学的環境によって変化する。動的電荷平衡法では、電荷に対し仮想的な質量M Qを導入し、電荷が次の古典的運動方程式に従うと仮定することによって、系の電子分極とその揺らぎを表現する。4）
これは、電荷が分子内で中性を保つという束縛条件を課したオイラー・ラグランジュ方程式から導かれる。もに関する和は、その分子内の原子サイトに関する和で、N nは分子内原子数である。個々の分子内ではなく、系全体で中性を保つとしても、同様の式が得られる。以上が通常の動的電荷平衡法の要約だが、このまま水中の電子移動系に適用すると、電荷の値が発散するという現象がしばしば起こる。これは次の単純な理由に依る。式(14)を用いた場合、全ポテンシャル(13)は電荷の組{Q}に関して２次式であるから、{Q}空間において負のエネルギー方向に発散する鞍型を取り得る。特に、水の様に正負に分極した原子間の静電引力が重要である系では、ゆらぎや衝突によって電荷が接近した場合に、このような現象が起こり得る。これは、式(13, 14)のモデルに本来的に内在する問題である。この問題を避けるためには、例えば、あるしきい値よりも電荷が増大した場合に跳ね返る様なポテンシャル壁を付け加えるということも考えられる。しかしそれでは電荷の軌跡に不自然な飛びが生じてしまう。そこで我々は、次の様なポテンシャル関数を提案した。51
これは、一見複雑に見えるが、余分なパラメータを付け足したりしていない点に着目して頂きたい。Qhに関して展開すれば、
となり、χ畠およびψは、元の物理的意味を基本的に保っている。ただし、後述のようにψには補正の余地を残している。この式(16)は水茶原子に関するもので、右辺第２項のために電荷QHは+2を超えることはない。また、全ての原子種に適用できるように関数型を一般化することも出来るが、実際問題としては、水の場合には水茶原子だけに式(16)を導入し、酸茶原子には元の２次式(14)を用いても、十分に安定なシミュレーションが実現される。同様に、水茶結合を形成するような極性の水茶原子だけにこの修正モデルを導入すれば十分である系は多いと予想している。式(13)中のクーロン項Jab(rab)には、簡単のために大野近似を用いた。勿論、オリジナルの崔荷平衡法の様に、この項を原子軌道関数間のクーロン積分として正確に計算しても良い。しかし、今回の研究においては、計算コストを上げてもそれをするべき積極的な理由は見出せなかった。
3.2 パラメータ決定の手続き
各原子における電気陰性度と硬さという、少数かつ物理的意味の明快なパラメータで系の電荷分布を決定出来ることが、電荷平衡法の魅力の一つである。今回は、酸茶原子のパラメータは文献6)から取り（χδ = 8.741 eV, ηδ = 6.682 eV)、水茶原子のパラメータは、気相中の水分子の双極子モーメントと分極率を再現するように最適化した。ここで、上述の沿の補正について触れておく。Emod (Q)では、Equad (Q)と異なり、Qに関する2次微分は定数ではない。仮にη 0 = η 0とすると、Q > 0における実効的な硬さ∂² E mod (Q) /∂Q²はη ⁰よりも大きくなる。 γ ⁰に補正の余地を残したのは、これを緩和するためである。ここでは、系の電荷分布を決定する主因子は1次微分すなわち電気陰性度であることから、ある特徴的なQの値において、Equad(Q)とEmod(Q)の1次微分が等しくなるという条件によって、ψを決定した。今の場合は、水分子中の水茶原子が問題になっているので、気相中の水分子の双極子モーメント（1.85デバイ）を再現する水茶原子上亀荷Q器（TIP3Pモデルの分子構造では、0.3287e）を用いて、QH = QHにおいて両者の1次微分が等しくなるとした。
3.3 分子シミュレーションによる再配置エネルギーの計算
第２節で議論した様に、水溶媒中の外圏型電子移動では、式(4)が成り立つ場合が多い。このとき、洛媒の電子分極は溶質の電子状態に断熱的に追随するので、分子シミュレーションにおいては、溶媒電荷の組を２組用意し、それぞれを電子移動の始状態と終状態に対応させることにする。それらを{Q t }、{Q f }とする。電子移動反応の分子シミュレーションにおいて鍵となる量は、垂直エネルギーギャップ
である。ここで、H iとH fは、反応の始伏態と終伏態における原子核運動を表すハミルトニアンで、それぞれのポテンシャルがV i (r, Q i )とV f (r, Q f )である。両者には異なった溶媒電荷の組が対応している点に注意して頂きたい。これら{Q'}、{Qf}は、与えられた原子核配置においてそれぞれ独立に、電子移動系の始状態および終状態の電荷分布に対して最適化されるものである。溶媒再配置エネルギー入は、次式から計算される。
ここで、〈…〉i,fは、始状態および終状態における熱統計平均を表す。図１は、シミュレーションから得られた再配置エネルギーを、浴質問距離R D Aの逆数に沿ってプロットしたものである。電子分極効果を含まない通常のTIP3Pモデルによる計算と、電荷平衡法によって電子分極を取り入れたモデルによる計算が比較されている。シミュレーションから得られた入は、誘電連続体モデルによる式(3)と同様、1/R DAに沿って直線によく乗っており、その傾きと絶対値は、電子分極効果を取り入れることによってスケールダウンされている。これは、定性的には式(3)の因子(1/ϵ∞ − 1/ϵo)で表現される現象である。式(2)からも明らかな様に、入は電子移動速度の活性化エネルギーを決定する因子であり、図に見られるような0.1 eVオーダーの違いも有意に効いて来る。ただし、次節に見る様に、量子効果を取り入れたFC因子の計算においては、入の値の違いを表面から消すような、あるスケーリング則が成り立つ。
3.4 久保-豊沢公式によるFranck-Condon因子の計算
式(1)の速度定数を時間領域の積分形式に変換し、キュミュラント展開を取ることにより次式が得られる。²）（ガウス型熱浴では高次のキュミュラントは消える。）
これは、高温極限で式(2)に帰着する。すなわち、△Vの平均とその揺らぎの時間相関関数を量子力学的に計算できれば、溶媒運動の量子効果を取り入れた電子移動速度が得られる。量子力学的な相関関数を具体的な分子モデルについて計算するのは一般には容易ではないが、ここではまず分子動力学シミュレーションから古典的な相関関数を計算し、そのフーリエ変換についてある量子古典対応関係を仮定することにより、それを求めるという方法をとった。7この対応関係は調和振動子系について厳密であり、溶媒中で熱平衡にある（すなわち強く励起されていない）結合モードの記述においても良い近似であると期待できる。調和近似が特に良くないのは遅い流体的運動であるが、それらの量子効果は大きくない。図２はシミュレーションから得られた古典的時間相関関数とそのフーリエ変換である。相関関数に見られる振動は800 cm−1周辺のピークに相当し、これは水溶媒の捻れ振動(libration)から来る。図に見られるように、電子分極を入れた電荷平衡法モデルではこの振動の強度が抑えられている。これは、溶質の電子状態変化へ応答する際に、固定電荷モデルでは原子核運動が全てを請け負わなければならないのが、分極モデルでは応答の一部を電子が引き受けるからであると解釈できる。図３は、反応速度定数の対数で定義した量子力学的活性化エネルギー △F*の△F°依存性である。古典極限である点線からのずれが量子効果を表す。いわゆる逆転領域△F 0 < −λにおいて量子効果が顕著であることが見て取れる。ここで、電荷平衡法と固定電荷モデルとの間で結果に差がないように見えるのは、図の軸に示されているようなスケーリングを行ったからであって、実際には図1に示した入の違いにより、活性化エネルギーすなわち反応速度の値は有意に異なる。見方を変えると、図３のスケーリング則が良く成り立つことを利用すれば、モデルを多少変更した場合を調べる際にも、式(20)以下に説明した計算を一々繰り返す必要はなく、再配置エネルギー入とX’だけを計算すれば、図３から量子効果を取り入れたketを直ちに見積もることが出来ることになる。なお、X′ ≡ 〈δΔV²〉/2kBTはΔVの熱揺らぎ幅から定義される再位置エネルギーで、調和振動子の場合はλと一致するのでしばしば混同されるが、本来は異なった量である。
4　おわりに
以上、竜荷平衡法を用いた分子動力学シミュレーションによって電子移動反応系における溶媒電子分極効果を調べた仕事の概要を紹介した。本稿では要点を絞るために詳細の多くを省略した上に、参考文献は最小限に留めてある。必要に応じて文献5.8）およびその中のリストから辿って頂ければ幸いである。近年発展の著しい第一原理分子動力学法91でも、原子数が千のオーダーの系を熱統計的に扱うのは、計算コスト的に未だ困難である。電荷平衡法の活躍する機会は多いと期待できる。第３節に述べたような方法論で他の色々な系を調べること10）や、そのためのパラメータの整備と改良は、今後も推進すべき課題である。ただし、対象とする系を拡げる際には、第２節に述べたような理論的基礎に注意を払いながら行うことが肝要であろう。さらに、式(4)の成り立たないような電子移動系のシミュレーションを行うにはどうしたら良いか、そこでは電荷平衡法を応用することが出来るのかどうかという問題は、断熱近似を超えた電子と核の動力学をどう記述し現実的なシミュレーションに載せるかという問題であり、これから理解を深めて行くべき課題である。