分子科学アーカイブス
AC0004

溶媒効果と MO 計算

公開日2007 年

佐藤啓文

著

6月

29 日

第１版＊

分子科学会編集委員会は、優れたテキストを分子科学アーカイブスとして公開しますが、その内容の一切の責任は著者にあります。読者からの貴重なご意見は、（edit-office@j-molsci.jp）で随時受け付けております。ご意見は編集委員会から著者にお伝えし、テキストの内容に反映していきます。
著者紹介佐藤啓文（さとうひろふみ）
所属：京都大学工学研究科分子工学専攻

分子理論化学講座

専門分野：理論化学・溶液化学・物理化学＊

第１版については、分子科学研究所・吉田紀生博士に精読して頂き、ご意見を頂戴いたしました。御礼申し上げます。
2007 年 6 月
分子科学会編集委員会

第四分科会溶媒効果と MO 計算

分子科学研究所佐藤啓文

hirofumi@ims

第四分科会・佐藤啓文

改訂版公開にあたって
1998 年の第 38 回分子科学若手の会夏の学校では、事前に参加者全員でメーリングリス
ト等を用いた討論を行い、内容について幾つかの誤りを指摘頂いた。本稿はこれらを反映し、実際の「夏の学校」時に配布された冊子体から若干の改訂を行ったものである。分科会は、当時広島大学の大学院生であった多田朋史博士（現：東京大学）および山崎健博士（現：アルバータ大学）が中心となり、メーリングリストは渕上壮太郎博士（現：横浜市立大学）がその運用をして下さった。また分科会の参加者全員が内容についての議論に参加して下さった。その名を付記することで謝辞としたい。
私自身が学生の時に「夏の学校」でよく耳にした「教科書以上論文未満」を意識して書いたのだが、時を経た今改めて見直すと、分かりにくい点や意図が伝わりにくい点、書き直すべき点も散見される。アーカイブスという性格を鑑み、あえてこれらに関しては訂正・加筆を行わなかったが、どうかご容赦願いたい。
2007 年 1 月
京都大学 工学研究科 分子工学専攻
佐藤啓文
第 38 回分子科学若手の会 夏の学校 第四分科会参加者
今井隆志（D2、分子研、平田研）
岡田一俊（D2、分子研、岩田研）
本田康（D2、京大、中辻研）
太田義啓（D1、北大、田中研）
岩橋建輔（D1、名大、大峰研）
帰家令果（M2、京大、梶本研）
向山綾子（D1、京大、梶本研）
秋元伸一（M1、筑波大、菊池研）
山崎健（M2、広大、高分子化学）
稲田久美（M1、広大、高分子化学）
山本亮介（M1、東北大、三上研）
光武亜代理（D2、分子研、岡本研）
篠田恵子（D2、東工大、岡崎研）
神坂英幸（D1、分子研、中村研）
豊田和男（D2、京大、中辻研）
関口健太郎（D1、京大、梶本研）
江頭和宏（D1、京大、梶本研）
渕上壮太郎（M2、東大、染田研）
杉木真一郎（M2、広大、高分子化学）
岸崇史（M1、東北大、三上研）
大森努（M1、京大、廣田研）
原野雄一（D2、分子研、平田研）
小林千草（D1、名大、大峰研）
多田朋史（D1、広大、高分子化学）
中田真秀（M2、京大、中辻研）
藤巻英司（M1、東北大、三上研）（順不同、所属と学年は当時のもの）

1

第四分科会・佐藤啓文

はじめにこの分科会では、理論化学の最近の手法、特に溶液内の化学過程を扱うための電子状態理論と分子論的な描像に基礎を置いた統計力学の手法について皆で勉強してしていこうと思う。今日、様々な化学過程に於ける分子レベルの微視的状態を知るための量子化学的アプローチや統計力学的アプローチは、それぞれの分野における発展を基に、統合され、より複雑で実践的な化学系へむかって発展を遂げている。これまでは複雑すぎて「理論化学」
では全く歯が立たなかった様なごくごく日常的な化学反応について、実際に観測されている多くの結果と理論計算の結果を直接つきあわせて議論をすることが可能となりつつある。
現時点での溶液内化学過程に対する理論は、過渡的要素を含みすぎていて普遍的な知識にはならないかもしれない。しかし、理論化学の目指すべき方向は、とりもなおさず化学の目指す方向の一つでもあるはずである。この分科会での勉強がそういったことへの考察の一つのきっかけになればと思う。
本稿は以下の様に構成してある。

1. 概観
2. 溶液状態を記述するための道具• 溶質・溶媒分子間の相互作用 — 分子軌道法とその周辺• 統計的手法による溶媒の分布決定 — RISM 理論
3. 溶液内分子の電子状態の記述 (1) 連続誘電体モデル
4. 溶液内分子の電子状態の記述 (2) QM/MM 法
5. 溶液内分子の電子状態の記述 (3) RISM-SCF/MCSCF 法
6. まとめと文献

最初に全体の概観について述べる。次に溶液の理論的記述の基礎となっている量子化学・統計力学の手法について簡単に説明をする。最後にこれらの手法を基に現在研究の進められている方法について概説する。

2

第 1 章 概観
1.1

溶液内化学過程を扱うための様々な手法

溶液内における化学過程を理論的に考察しようとした場合、二つの観点が重要になってくる。一つは量子化学の観点である。多くの場合、化学結合が新しく出来る／壊れる、注
目する分子が光励起をする、といった事象が化学的な興味の対象である。これらの現象に伴って変化する電子状態の記述（主に分子軌道法に基づく）は、殊に孤立分子の場合は、今日の理論化学における中心的課題の一つである。二つ目は、溶液全体の自由度の観点である。1023 個もの分子の相互作用の結果として、分子集団としての溶液の性質が発現してくるわけであるから、そのままの溶液全体を考察するために統計的手法が必然的に重要になってくる。このような取り扱いの多くは古典力学に基づいているが、物理化学の問題の一つとして古くから沢山の研究がなされてきている。昨今は計算機の発達により、化学過程を計算機上で数値的に実験してみるという分子シミュレーションが盛んである。
もちろん、注目する溶質分子と周辺の分子との相互作用が極めて小さい場合は、溶質分子を孤立した分子と見做して考察することがよい近似になるだろう。また、系の（化学的に）静的な熱力学的性質を問題としているときは、逆に統計力学的考察で十分であることが多い。残念ながら（かどうかは意見が分かれるかもしれないが）、多くの化学過程はその双方を考慮しないとうまく記述できないようである。そこで、より実戦的理論として両者を組み合わせたアプローチがいくつか提案されている。
Table 1.1 に本分科会の内容に関連する代表的な手法の幾つかをまとめた。おおよそではあるが、縦のならびは量子論的側面に、横のならびは自由度の取り扱い方に対応している。
表 1.1: 溶液内化学反応を調べるための手法

classical model↑

連続体モデルPoisson eq.
Born model

SCRF, MPE,
↓ IC, COSMO,
quantum model PCM, GB, etc

シミュレーション
MD (mol. dynamics)
MC (Monte Carlo )液体の統計力学
MSA
RISM

QM/MM

RISM-SCF/MCSCF
略称に関しては後述する。

まず、一つ目の観点、量子力学的側面を考える。最も単純な考え方は、なんらかの方法

3

第 1. 概観

第四分科会・佐藤啓文で決定された分子の電荷分布がつくる古典的電場のみとすることである（つまり量子論的考慮は全くしない）。量子論的考察をするためには、これまで蓄積されてきた孤立分子系
での手法を取り入れるのが適当であろうから、半経験的手法や非経験的な手法が用いられている。
二つ目の観点、自由度に対してはどうだろうか。これまでに提案されてきたアプローチ
は主に三つの立場に分けることができる。一つ目は、「無数の溶媒分子を個々の分子として扱うことをせずに、なんらかの方法でより簡単化する」立場である1 。なかでも電子状態理論との組み合わせが最も盛んなのは連続体近似である。この方法は溶媒を大胆にモデル化して個々の分子構造は無視してしまい、連続体として取り扱ったものである。通常は静電的相互作用が重要であることが多いので、誘電体として扱う。後に述べる Generalized
Born 法の様に手軽に計算が行える半経験的分子軌道計算で実験結果を比較的よく再現できることもあり、現在でも最も広く用いられている。しかし溶媒分子の個性の全てを誘電率で表現してしまおうという立場を取るこの方法では明らかに不十分であることが多く、より現実的な溶媒モデルの模索が行われている。その一つが、二つ目の立場である、「統計力学的側面をシミュレーション（計算機実験）に基づいて考慮する方法」である。三つ目の液体論は「同じ統計力学的側面を、計算機実験ではなく、純粋に理論的に考察しようとする」立場である。いずれの方法も、量子化学計算との結合がなされていて現在種々の化学系への応用がなされている。本分科会ではそれらの手法について見ていこうと思う。

1.2

平均場近似と単純化したモデル

本題に入る前に、溶液系全体の取り扱い方についてもう少し補足しておく。平衡状態にある無限希釈系の溶質分子を考えているときに、最終的に我々が必要なのは注目している
溶質分子の波動関数である。波動関数が得られれば、エネルギーはもちろんのこと分子の電子状態や周辺との相互作用を理解することができるだろう。Born-Oppenheimer 近似の
もと、溶質分子の核座標を R として固定し、電子の座標を r とすると、孤立した状態にある分子のエネルギー（Evac ）は電子ハミルトニアン（H）を使って、∫Evac (R) =drΨ∗0 (R, r)H(R, r)Ψ0 (R, r)= < Ψ0 |H|Ψ0 >,(1.1)と書ける。実際には系すべての分子を量子論的に扱うことができないので、溶媒部分は古典的な取り扱いをせざるをえない。従って与えられた配置に対しての全体のエネルギーは
溶質・溶媒間の相互作用 V と溶媒の寄与 Esolv を使って、Etotal =< Ψ|H + V|Ψ > +Esolv ({R0 }),(1.2)となる。このときの V は、溶媒分子の核座標の組 {R0 } と溶質分子の核座標 R、波動関数|Ψ > で評価される。

V ≡ V({R0 }, R, Ψ),

(1.3)

1

例えば周辺の数分子のみを考慮して溶液のモデルとする「超分子法」は、近接する分子の自由度のみを考慮する、という点ではこのカテゴリーに属すると考えてもよいだろう。しかしこれは溶液というよりはクラスターのモデルである4

第 1. 概観

第四分科会・佐藤啓文

系全体の自由エネルギー Gsol は、溶媒の配置のアンサンブルをとる（< ... >Γ で表す）ことで決まる。

1
Gsol = − ln < exp[−βEtotal ] >Γ .
β

(1.4)

注意しなければならないのは、これらの式に含まれている |Ψ >、{Γ} は相互に関連していることである。つまり溶質分子の電子状態が変化すれば、それの作る場は変わるので溶媒の構造（配置）も変化する。逆に溶媒の構造が変化すれば、溶質分子上にできている場も変わるので、その電子状態も変化するわけである。つまり平衡状態では Ψ、{Γ} を含め全体が無矛盾 (self - consistent) である必要がある。
平衡状態を考えているときなどは、式 (1.4) で溶質の核座標は固定されていているので、

1
Gsol = − ln < exp[−β(< Ψ|H|Ψ >)] exp[−β(< Ψ|V|Ψ > +Esolv )] >Γ
β1= < Ψ|H|Ψ > − ln < exp[−β(< Ψ|V|Ψ > +Esolv )] >Γ
β

(1.5)

つまり、溶媒全体が溶質分子に及ぼす影響をその平均で考えることができる。これが平均場近似（mean ﬁeld approximation）である。この場合でも |Ψ > と {Γ} は相互に関係しているので self - consistent な取り扱いが必要である。溶媒の作る場は時間スケールで見ると、それらの分極的効果による非常に応答の早い寄与と、分子そのものの運動（回転など）に起因する応答の遅い部分からなっており、平均場近似ではこれらがいずれも緩和していることを要求している。従って、溶質分子の電子状態の変化に対して溶媒が遅れてゆっくりと応答する様な系を扱うことは一般には難しいということになる。
計算時間を節約するための近似として、あらかじめ相互作用 V を波動関数に依存しないとして与えてしまう方法が考えられる。例えば、量子化学計算を孤立している溶質分子についてあらかじめ行って、これから相互作用 V を決めてしまえばよい。

Gsol ≈ Gapprox
1= − ln < exp[−β(< Ψ|H|Ψ > +V + Esolv )] >Γ0
β1= < Ψ|H|Ψ > − ln < exp[−β(V + Esolv )] >Γ0 ,
β

(1.6)

このとき、溶媒のアンサンブルも {Γ0 } になってしまうことに注意してほしい。勿論、この方法では |Ψ > と {Γ} について self -consistent な取り扱いにはなっていない。
Blair らは水中のホルムアルデヒドの吸収スペクトルのシフトを調べるために、この方法を用いた [12]。彼らはまず、真空中の量子化学計算を行って溶質分子（ホルムアルデヒド）の電子状態とそこから導かれる有効電荷を決定した。得られた結果から相互作用 V を決め、通常の MD 計算を行って溶媒分子の配置のアンサンブルをとる。さらに、これらの分布の中から幾つかを選んでその配置の点電荷の作る静電場の下で、再度量子化学計算を行ってスペクトルシフトの理由を説明した。
この方法は、計算時間を大幅に削減できる点において非常に有効で数多くの例があるのだが [13]-[15]、すでに指摘したように全体で無矛盾にならない。上式中の V の |Ψ > によ5

第 1. 概観

第四分科会・佐藤啓文

る評価の仕方は方法によって微妙に違うが、いずれも孤立分子（溶質）の波動関数を用いる。最終段階でスペクトルを計算するためには溶媒からの寄与を含んだ溶質分子の波動関数を求めているのだが、溶媒構造（アンサンブル）は従来の古典論のシミュレーションで得られる結果の域を出ないし、溶質の波動関数も適切に改善されているわけではない。
この分科会の論点である現在の溶液内電子状態理論は、平衡・非平衡状態における溶質
分子の電子状態を求める方法であるので、孤立分子系の電子状態に基づいているこれらの研究については、これ以上言及しないことにする。

6

第 2 章 溶液状態を記述するための道具
前述したように、今日の大部分の溶液の電子状態理論が対象としているのは無限希釈系のみである。しかも溶質分子は系の中に唯一つ、通常は座標固定された分子として存在しているとする取り扱いが一般的である。溶媒分子は周辺に多数存在していて、溶質分子の電子状態にある種の外場として作用している。この際、溶媒分子の核の運動については古典的に取り扱うこととする。溶質分子の核の運動は、他の全ての自由度に対して断熱的記述が適当であると、ひとまず仮定しておく。
ここまで問題を限定しても、溶媒分子は事実上無限数存在し、しかもそれらの配置に対しての熱平均を取る必要がある。従って、現段階の理論を構築していく上での問題点は、それぞれの溶媒分子（溶質分子・溶媒分子間の相互作用）の記述を如何に単純化するか、
そしてそれらのアンサンブルをどう取るか、の二点に集約されるだろう。
このうち、前者は溶質・溶媒間の相互作用ポテンシャルの記述の問題と言い換えることもできる。また後者については、本質的に統計的考慮が必要であり、シミュレーションあるいは積分方程式といった手法を使う（最も広く使われている連続体近似という考え方は、
これら溶媒構造を極端に単純化あるいは無視したことに相当すると言える）。この章では、これらの観点から、これまで提案された溶液内分子の電子状態を調べていくための手法について解説していく。

2.1

溶質・溶媒分子間の相互作用ポテンシャル

2.1.1

溶媒分子のモデル化

これまで提案されてきた、溶質・溶媒分子間の相互作用の記述について、代表的な研究の幾つかを大凡に分類したのが Table 2.1 である。
Table にもあるように、これまでの手法の中で完全に量子化学的な取り扱いをするのは、超分子法だけである。この方法では、注目する溶質分子と周辺の幾つかの溶媒分子とから成るクラスターを考えて、全体を標準的な量子化学計算で扱う。こうした方法が溶液モデルとして用いられることが多かった。勿論、こういったモデルでも溶媒分子の個数を十分増やして、シミュレーションの手法等でそれらの作るポテンシャル面をつぶさに調べ上げていくことが出来れば良いのだが、多数粒子系に対しては現実に計算を行うことは不可能である。限られた個数の溶媒分子のみを考えることで計算は可能になるのだが、本来無数の溶媒から働く遠距離力（クーロン力）や、溶媒分子のアンサンブルを取るうえでの困難など、あくまでクラスターの計算であって、溶液系とは区別されるべきものである。

7

第 2. 溶液状態を記述するための道具

第四分科会・佐藤啓文

表 2.1: 分子間相互作用の記述

full quantum↑

超分子法 ≡ 量子化学的取り扱い
擬ポテンシャル法

Honda-Kitaura ポテンシャル↓
fulll classical12-6-1 型ポテンシャル, 多重極子展開

こうした全系の量子力学的な取り扱いが出来ない以上、溶質分子の電子状態についての
み量子化学的に取り扱い、溶媒を古典的に取り扱う、いわゆる混成的な方法が現実的であり、主流である。
なお、本節に挙げたモデルの全てが必ずしも溶液内分子の電子状態理論を目的として提案されたものではないことを、予め注意しておく。

擬ポテンシャル法
量子化学的手法との組み合わせの中で、溶媒分子の簡単化に際して最も直接的なのは擬ポ
テンシャル法であろう。現在広く用いられている、重原子に対しての eﬀective core potential
法・model potential 法 [18] 等では、化学的特性に直接関与が小さいと考えられる内殻電子については、価電子に対しての eﬀective なポテンシャルであるとして、電子状態を解く際に予め与えてしまう。この結果、計算量を削減することができ、非常に大きな成功を納めてきた。同じように溶質分子の電子状態を決定する上で、溶媒分子からの直接的寄与が一種の外場として働く eﬀective ポテンシャルに置き換えることができれば、計算時間を大幅に削減できる。これにより簡単に溶媒分子の総数を増やすことが可能になり、シミュレー
ションの手法と組み合わせることができるだろう。
model potential 法の拡張として分子系への最初の応用がなされたのは、Ohta らによる
eﬀecitve fragment potential 法 [19] である。彼らは NH3 を二電子系として扱い、残りの電子は擬ポテンシャルに置き換えた。実際に幾つかの分子間でのポテンシャルを報告しているが、全電子で解いた時に比較して非常に良い結果を得ている。
直接的に溶液内電子状態への拡張を目的とはしていないが、Katsuki らによるスペクト
ル表示の方法 [20] も、大規模系の記述に際して非常に強力であると思われる。この方法では Fock 演算子 F を次式のように active な部分と frozen な部分に分割する；


∑ ZA

1
F = − ∇2 + −
+
2
A rA


f rozen∑+ Ωactive∑(2Jj − Kj ) − Vad j

(2Jj − Kj ) + Vad  Ω,

j

8

(2.1)

第 2. 溶液状態を記述するための道具

第四分科会・佐藤啓文

ここで、Z は核電荷、J 、K はクーロン、交換積分で Vad は適当な演算子である。Ω は任意の演算子 O に対して、基底 {χa } の重なり積分 S を使い、

ΩOΩ =

∑

|χa > (S −1 )ab < χb |O|χc > (S −1 )cd < χd |,

(2.2)

abcd

と定義され、これを演算子 O に対してのスペクトル表示と呼ぶ。もし、この基底が完全系であれば上の Fock 演算子は厳密に正しい。彼らはこれをアンモニアの系等に適用して、良い結果を得ている。
本格的な溶液系を志向した擬ポテンシャル法を用いた最初の例は、Vaidehi らによってなされた研究であろう [21]。彼らは擬ポテンシャルで表された溶媒分子モデルを用いて、abinitio 法と MD を組み合わせた。彼らの Fock 行列要素は

Fµν = Hµν + U

PP

∑[

]

1
(µν|λσ) − (µλ|νσ) ,
+
2
λ、σ

(2.3)

である。ここで、 U P P が各溶媒分子から来る擬ポテンシャルで、

U

PP

(

∑

)

AB
qB
1
(ri ) =++ AB−
exp(−αB |ri − rB |2 ),
2
|r
−
r
|
|r
−
r
|
i
B
i
B
B

(2.4)

を用いる。式中のパラメーターは全電子の計算結果を再現するように決定する。ポテンシャ
ル関数の形やパラメーターの決定等の問題点が残っているものの、最初に溶液系に取り組んだ意義は大きいだろう。彼らは最初 Li+ –H2 O 系の計算を行ったが、後に密度汎関数理論(DFT; Density Functional Theory) と結合して H2 O–H2 O 系などの研究も行っている。

Honda-Kitaura 型ポテンシャル
こうした溶質・溶媒間の相互作用の単純化をさらに進めた一例が、Honda-Kitaura のポ
テンシャルである [22]。彼らは、波動関数に依存する項をそれらの重なり積分の自乗のみで書き直してしまい、残りの相互作用を有効電荷（q ）間のクーロン力だけで表した。
、b
qr qs a∑
+
Vab =
Cij
r
i、j
r 、s rs

a
、b
∑

[∫

]2

φi φj dτ.

(2.5)

ここで C はパラメーター、φ は局在化した分子軌道である。形式的に Murrel の double
perturbation approach[23] に等価と思われるが、上述のスペクトル表示法とも関連づけられるだろう。この方法は電子相関を含んだ場合など、非常に多くの系に応用されており、
Muguruma らにより遷移金属錯体への応用、三体以上の多体表現を含む方法についても研究されている [24]。

9第 2. 溶液状態を記述するための道具

第四分科会・佐藤啓文

古典論的ポテンシャル
さらに単純化を押し進めるには、Honda-Kitaura 型ポテンシャルの重なりの項を、なんらかの（一般的には距離に依存する）解析関数に置き換えてしまえばよい。これは、もはや量子化学的描像は暗に含まれているとした古典的相互作用ポテンシャル [25] であり、MDや MC のシミュレーションや積分方程式論で広く用いられる。代表的な例は Lennard-Jones型の関数で表現する 12-6-1 型ポテンシャルであり、分子 A と B の相互作用は、[(∑

EA−B =

4²a,b

a⊂A、b⊂B

σab
r

)12

(

σab
−
r

)6 ]

+

qa qb
r

(2.6)

とかかれる。この表式では、分子間の相互作用を幾つかの代表的な点（a、b : 原子などをとる）の間の相互作用の和で近似する。ここではこれ以上触れないが、12-6-1 型の相互作用関数が、全ての化学系に対してよい近似になっているわけではないことを注意しておく。

2.1.2

分子のモデル化と静電ポテンシャル

さて、ここで少し見方を変えて、分子の電荷分布 ρ(r) とそれがつくるポテンシャル場について考えてみる。周辺に存在している分子との相互作用のなかで静電的寄与は中心的な役割を果たすからである。元の量子化学的な考えかたに沿えば、電子波動関数 Ψ(r) のつく∑
る電荷分布は密度演算子 ρ̂(r) = N
a=1 δ(ra − r) を用いて、

ρQM (r) = − < Ψ|ρ̂(r)|Ψ > +N uclear∑

ZA δ(rA − r),

(2.7)

A

と書くことができる。12-6-1 型などの点電荷を用いる古典論的ポテンシャル関数の場合は、

ρP C (r) =site∑qi δ(ri − r),

(2.8)

i=1

となる。ここで導入した有効電荷（qi ）は何らかの手続きで決めなければならない。これ
らの電荷分布が点 r0 に作る静電ポテンシャル V (r0 ) は、

VQM (r0 ) =
0

∫

dr

ρQM (r),|r − r0 |

dr

site
qi
ρP C (r) ∑=,0|r − r | i=1 |ri − r0 |∫VP C (r ) =(2.9)と書ける。同じことを多重極子展開を基に書くこともできる。Laplace 方程式の一般的な解は極座標表示を用いるとVM P E (r, θ, φ) =

∞ ∑
+n (
∑
n=0 m=−n)Bnm
Anm r + n+1 Pnm (cos θ)eimφ ,rn(2.10)となる。ここで Pnm (cos θ) は Legendre 陪関数であり、Anm は定数 Bnm は電荷分布に関連する定数である。

10

第 2. 溶液状態を記述するための道具

第四分科会・佐藤啓文

8

8

(b)

(a)

-8

-8
8

-8

8

-8

8

8
(d)

(c)

-8

-8
8

-8

-8

8

図 2.1: 水の分子モデルのつくる静電ポテンシャル
それぞれ、(a) は双極子モーメント（µ=2.2D）、(b)TIP3P モデル（12-6-1 型ポテンシャル関数）、(c) 非経験的分子軌道法（Hartree - Fock (HF) 法、STO-3G 基底関数）、(d) 非経験的分子軌道法（Hartree - Fock

(HF) 法、6-31G**基底関数）により計算したもの。座標軸は au 単位で酸素原子を原点とし、等高線は-0.5
から+0.1 まで 0.01au/e 毎である。

水分子のモデルについて、これらのポテンシャルを表したのが Figure 2.1 である。
図のいずれも、中心の太い縦線左側がポテンシャルが正の領域、右側が負の領域になっている。分子からある程度離れた領域ではどれも極めてよく似た結果を与えているが、近い領域では大きく異なる。二つの分子軌道法をもとにした結果を比較しても、より大きな11

第 2. 溶液状態を記述するための道具

第四分科会・佐藤啓文

基底関数である (d) は酸素付近のポテンシャルが負の領域が浅くて広い。但し、この図は静電的な寄与しか考慮していないことに注意してほしい。分子近傍ではそれ以外の効果（交換相互作用、電荷移動など）も重要であり、実際の分子間の相互作用はもっと複雑である。

2.1.3

分子軌道法を通した分子の見方

分子のモデル化という観点で、分子軌道計算により得られた結果は含まれている近似の正体もはっきりしているので、理論的信頼がおけそうである。しかしながら、ひとくちに非経験的な分子軌道法と言ってもその近似の精度は様々であり、適切な手法を用いないと定量的には勿論、定性的にも誤った結果を与える。よく知られているキーワードは「電子相関」と「基底関数」であろう。当然のことであるが、両者ともより高い精度まで取り扱えることが望ましいが、現実には計算時間の制約があって思うようにはいかないだろう。
本節では分子軌道法から得られる結果について、実際の計算をもとに議論する1 。
例として、CO の双極子モーメント (µ)（Table 2.2）と水の平衡構造 (re 、θe )・双極子モーメント (µ)（Table 2.3）の計算結果をまず見てもらう。

method

表 2.2: CO の双極子モーメントµ/a.u. method

HF / STO-3Ga
HF / 4-31Ga
HF / 6-31G*a
HF 極限 a
実験値

0.066 153 次元 CIb-0.237 2484 次元 CIb-0.131 11 次元 MCSCFb-0.1100.044µ/a.u.
0.048
0.130
0.066

A.Szabo and N. Ostlund、「新しい量子化学」から抜粋
aA. D. McLean and M. Yoshimine,

Int. J. Quantum. Chem. 15,313(1967)b

岩田末廣、分光研究、30, 3(1981) から抜粋

CO の双極子モーメントは、正しい実験結果によれば炭素の側が負極である。電気陰性度の考えに基づいた直感では酸素側が負になりそうだが、この場合その他にも重要な寄与がある。炭素の lone pair が結合とは逆の方向を向いていることである。二つの逆方向に働く寄与は拮抗しているのだが、後者が僅かに勝るので実験で示される様な結果になる。表に示されているように、STO-3G 基底以外の Hartree - Fock (HF) 法は軒並み誤った符号を与えている。しかし、電子相関をとりこむことによって（CI/MCSCF 法）、符号は改善さ
れている。

1

分子軌道計算に慣れ親しんでいない人には少々辛い内容かもしれないが、分科会でじっくり取り組めばよい

12

第 2. 溶液状態を記述するための道具

第四分科会・佐藤啓文

表 2.3: 水分子の構造と双極子モーメント

re /Å

θe / deg. µ/ a.u.
re /Å θe / deg. µ/ a.u.
Hartree-Fock2nd. Møller-Plesset
STO-3G 0.9901
100.0
0.679 1.0139
97.2
0.652
4-31G
0.9509111.2
1.026 0.9747
108.8
0.994
6-31G*
0.9478105.5
0.876 0.9689
104.0
0.869
6-31G** 0.9430
106.0
0.860 0.9610
103.9
0.825
re /Å θe / deg. µ/ a.u.
39-STO / SCF a
0.9398
106.1
0.785
実験値

re =0.9573 Å

re /Å θe / deg. µ/ a.u.
39-STO / SDCI a
0.9525104.9
0.755

θe =104.5 deg.

µ= 0.728 a.u.

A.Szabo and N. Ostlund、「新しい量子化学」から抜粋
aB.J.Rosenberg and I. Shavitt, J.Chem.Phys., 63,2162(1975)

次に水の例に移る。Table 2.3 に示されている、Rosenberg らによる結果は 39 個の Slater型拡張基底関数を用いている。この基底関数は HF 極限に近いエネルギーを与えている。
また SDCI については完全 CI による相関エネルギーの約 92%を与えていると推定されている。計算結果をつぶさに見てみると、平衡核間距離は HF 法では概して小さく見積もられていることが分かる。角度については二次の摂動論 (MP2) の計算はおおむね HF の誤差を小さくしている。双極子モーメントを見てみる。STO-3G は（少し小さめではあるが）、よい結果を与えているようだが、少し基底関数を大きくすると逆に大きくなる。39-STO が
HF 法で最良の結果だとすると基底関数を大きくするに従って、（STO-3G を除いて）双極子モーメントは徐々に小さくなっていくようである。電子相関を考慮することでも全体に小さくなっている。
さて、どう思うだろうか？
（今更、「やっぱり HF 法はダメで電子相関を入れなければ」と言うありきたりな答えは期待してない）
二つの点を指摘をしておく。第一に、目的に対して適切なモデルを組み立てていなければ計算としての意味はない、ということである。特に極端に精度の低い計算の場合、一見すると実験値と比較して適切な物理量を導き出しているように思えることもある（CO のSTO-3G の例など）。しかし様々な角度から計算として信頼できるものかを吟味する必要がある。これは電子相関を入れる場合の手法の選び方についても同じである。どんなに結果の数値がもっともらしく見えたとしても、例えば一重項ジラジカル的な電子構造を持つ解
離極限を MP2 で計算するわけにはいかないし、不適切な空間を選んだ CASSCF の計算は不適切な結果を与えるだろう。大切なことは数値があっているかどうかではなく、適切な計算モデルを組み立てているかどうかである。数値の合致は、その次の段階の問題である。
もう一点は、計算の評価の仕方である。例えば、「電子相関」の寄与の大きさを知るため

13

第 2. 溶液状態を記述するための道具

第四分科会・佐藤啓文には本来は用いた基底関数系における正確な値（つまり完全 CI の値）と比較すべきであって、実験値と比較すべきではない、ということである。ただ、残念ながら、今用いているすべての基底関数系に対しての正確な値は知られていないので、標準的指標（便宜）として実験値を用いざるを得ないのである。同じことは基底関数についても言える。基底関数
が十分であるかどうかは、その基底関数系を用いて得た物理量が HF 極限（あるいはそれに近い）を与える基底関数系を用いて得た物理量と同じかどうかで判断されるべきある。
CO の例でいえば、実験を再現するためには電子相関が本質的に重要であり、いくら基底関数を大きくしても HF 法による双極子モーメントの符号は誤った結果を与える。HF 法で（おなじ程度に）誤った結果を与える基底関数系は、適切な大きさだったからこそ誤った、ともいえる。
勿論、実験結果との比較に意味がない、ということを言っているのではない。実験結果との比較に入る以前の問題として、理論計算として適切であるかを十二分に検証する必要がある、ということである。

最後に、本分科会で最も重要で基礎的な分子間相互作用エネルギーについて BSSE(Basisset Superposition Error) という観点で考えてみる [27]。普通、二つの分子の相互作用エネ
ルギー ∆EAB を計算をするためには、二つの分子からなる系のエネルギーを計算し、それぞれの分子のエネルギーを引いて相互作用エネルギーとする2 。

∆EAB = EAB (ξA ⊕ ξB ) − EA (ξA ) − EB (ξB ).

(2.11)

ここで ( ) の中は計算で使う基底関数系を表している。個々の分子を計算する場合に比べて、二分子を計算しているときの方が関数系が張っている空間は広い。結局、二分子の計算の方が記述がよいことになるので、この方法で得られた相互作用エネルギーは安定化を過大評価しているはずである。当然のことながら個々の分子に対して用意した基底関数
が十分大きければ、この過大評価は無視できるほど小さくなる。
この効果を補正するためには、個々の分子を計算する際ももう一つの分子があるべき所に空の基底関数を用意して置けばよい。これが Counterpoise(CP) 法である [26]。
CP∆EAB= EAB (ξA ⊕ ξB ) − EA (ξA ⊕ ξB ) − EB (ξA ⊕ ξB ).

(2.12)

Figure 2.2 に Hartree-Fock 法による水の二量体のポテンシャル面を示した。CP 法の補正を行った結果も一緒に示してある。安定化エネルギーは約 5kcal/mol と見積もられている。STO-3G も一見するとよい結果を与えているのだが CP 法を行うと安定化エネルギーは約 2kcal/mol、ということになってしまう。一方、比較的大きな基底関数である 6-31G**
では CP 法による補正を行っても結果はそれほど変わらない。
分子間の相互作用は、(1) 静電的相互作用、(2) 交換相互作用、(3) 電荷移動、(4) 分極効
果、などからなっていると考えることができる。CP 法による補正は (3) を取り除いてしまうので、この寄与が大きな系で使うには不適切だが、今のように基底関数依存性が大きな
場合は、基底の張っている空間の広さの違いから来る、計算上の artiﬁcial な問題である。
2

ここでは Hartree-Fock 法のみを考えるので size-consistency とはまた別の問題である。

14

第 2. 溶液状態を記述するための道具

O

4.0

Relative Energy / kcal mol -1

(a)

H
H
H

O

4.0

H
2.0

∆E (without correction)∆E (with CP correction)
0.0

-2.0

2.0

∆E (without correction)∆E (with CP correction)
0.0

-2.0

-4.0

-4.0

-6.0

-6.0

2.5

3

3.5

(b)

6.0

Relative Energy / kcal mol -1

6.0

第四分科会・佐藤啓文

4

4.5

5

2.5

Distance O-O /Å

3

3.5

4

4.5

5

Distance O-O /Å

図 2.2: 水の相互作用ポテンシャルと BSSE(a) STO-3G による結果 (b) 6-31G**による結果
もとの安定化エネルギー（∆EAB ）という側面では二つの計算結果は約 5kcal/mol とそれほど変わらないが、その相互作用は全く異なった内容を持っていたことが窺える。つまりは、この計算に於ける STO-3G の結果は正しい物理的描像を反映できていないようである。
ごく当たり前のことだが、数値をいつでも鵜呑みにすべきではない。最も大事な点は、どんな物理量に注目しているのか、それはどんなモデル（計算手法）で得たものならば適切であるか、ということである。分子軌道計算の手法はそれぞれ明解な論理構造をもつ近似から成り立っており、注目する物理量を知るために必要な計算手法は自ずから決まってくるはずである。数値を再現することは勿論大事だが、それ以前に、論理的に・物理的に正しい描像でモデルを組み立てなければ理論としての存在意義がなくなってしまう。

2.2

溶媒の分布決定

溶媒の分布を決定するための統計的な手法としては二つの方法が考えられる。一つは、MD あるいは MC といったシミュレーションであり、今一つは相互作用点モデル ( RISM;
Reference Interaction Site Model) 等の積分方程式である。シミュレーションについては改めて説明する必要もないだろう。計算機を用いて溶液状態をいわば数値的な実験として再現する方法である。今日の溶液研究に於ては中心的な役割を果たしている。これまで述べてきたようなポテンシャル関数を選んで分子のモデルとし、MD 法では Newton 方程式に従って運動する分子の軌跡を時間に沿って追跡していく。この方法の最大の特徴は時間依存する物理量を直接観測できることにある。一方の MC 法は配置空間の中で乱数を発生し、平均値として物理量を得る手法である。いずれの方法も実際の計算では、これまで開発された様々な技術を基にしているが、直感的には非常に理解しやすいことと思う。

15

第 2. 溶液状態を記述するための道具

第四分科会・佐藤啓文

後者の積分方程式は、分子性液体に関する統計理論にその基礎をおいている。液体論は長い歴史と多くの輝かしい成功をこれまで納めてきている。液体論のなかでも「分子」の理論という点で、RISM は特筆すべきである。溶媒構造は溶質分子上に決定された各相互作用点からの一次元の動径分布関数に射影されるが、計算時間・得られた結果ともに非常に実際的である。とくに今の目的である、電子状態との関連を考える上では、誘電体近似を超える有力な手法として期待される。本節では溶液の構造を調べるための基礎的な統計力学の手法について平田の講義ノ−トを基に [9]、概説することにする。分科会の主な目的
は液体論の勉強ではないので基礎的な事項にとどめ、その詳細については、多くの他の著書を参照されたい [8]。

2.2.1分配関数と分布関数

この節では、まず多粒子系の記述の仕方、すなわち分布関数について学ぶ。分布関数は溶液の構造を表現していると同時に、多くの熱力学量を統計力学から導く上で極めて重要である。
N 個の粒子からなる系を考える。各座標を rN ≡ {r1 , r2 , · · · , rN } とすると、ハミルトニ
アン H は相互作用ポテンシャル U を使って、

H=

N
∑
i

1 2
pi + U(rN ),2mi(2.13)となる。カノニカル統計集団の分配関数は、β = 1/kB T を用いて、

QN =

∫
∫ ∞
∫
1 ∫∞
dr1 dr2 · · · drN e−βH ,
dp
dp
·
·
·
dp
·
·
·
·
·
·
1
2
N
N !h3N −∞V−∞
V(2.14)となる。このうち、運動量に関する因子はガウス関数であり、先に積分できる。

QN =

1
ZN .
N !Λ3N(2.15)

ここでドブロイ熱波長Λ=h√,
2πmkB T∫

配置積分 ZN =

≡

∫

V
V

···
···

∫

∫V

dr1 dr2 · · · drN e−βU(r1 、r2 、r3 、···、rN )drN e−βU(r ) ,N(2.16)Vである。系の Helmholtz の自由エネルギー A は、分配関数 QN を使って、

1
A = − ln QN
β(2.17)と表される。
液体中の分子は熱運動によってその位置を変えているが、ある瞬間における全粒子の位置を qN ≡ {q1 , q2 , q3 , · · · , qN } と書き、系の配置（conﬁguration）と呼ぶ。式 (2.16) の被16

第 2. 溶液状態を記述するための道具

第四分科会・佐藤啓文

積分関数はあるひとつの配置が出現する確率密度 PN (qN ) に比例しているので、配置積分
ZN は確率の規格化定数と見做すことができる。

e−βU(q )PN (q ) =
.
ZN
N

N

(2.18)

したがって、配置に依存する物理量を A(qN ) とすれば、そのカノニカルアンサンブル平均
は全ての配置 {qN } についての積分を取ればよいので、∫···< A >≡V∫

V

A(qN )PN (qN )dqN ,(2.19)と定義することができる。
ところで、空間のある位置における局所的な分子の密度（密度場）を定義するために、デルタ関数を用いることもできる。ある瞬間に i 番目の粒子が空間のある位置 ri にあったとする。この時、空間の任意の点 r の周りの微小体積要素 dr に粒子のある確率はδ(r − ri )dr,

(2.20)と表すことができる。同様に、一番目の粒子が r0 の周りに、二番目の粒子が r00 の周りに、
・
・
・
(N )
N 番目の粒子が r の周りにある確率は、

δ(r0 − r1 )δ(r00 − r2 ) · · · δ(r(N ) − rN )dr0 dr00 · · · dr(N ) ,(2.21)

で与えられる。この関数のそれぞれの粒子の位置 rN ≡ {r1 , r2 , · · · , rN } についての熱平均
値をとることで上述の確率密度 PN を定義することができる。

PN (r0 , r00 , r(3) , · · · , r(N ) ) =< δ(r0 − r1 )δ(r00 − r2 ) · · · δ(r(N ) − rN ) > .

(2.22)

従って熱平衡状態において、点 r0 のまわりに 1 番目の粒子を、点 r00 のまわりに 2 番目の粒子を・・・点 r(N ) のまわりに N 番目の粒子を見いだす確率は、

PN (r0 , r00 , r(3) , · · · , r(N ) )dr0 dr00 · · · dr(N ) ,(2.23)

で与えられる。
この分布関数 PN は分子の位置に関する全ての情報を含んでおり、原理的にはこれを基に全ての熱力学量を求めることができるはずである。しかし、このような関数を解析的に求めることは困難であり、より少ない自由度を含む計算可能な分布関数を求めることが実際的になってくる。そこで、N 個の粒子のうち、特定の n 個の配置 (r0 , r00 , · · · , r(n) ) に注目して、残りの (N − n) 個についてはどこにあってもよいと考える（つまり、この自由度を積分してしまう）。
(n)
PN (r0 , r00 , · · · , r(n) ) =

∫

V

···

∫

V

PN (r0 , r00 , · · · , r(N ) )dr(n+1) dr(n+2) · · · dr(N ) .

(2.24)

この分布関数は特性的確率分布 “speciﬁc (n-particle) distribution function”と呼ばれている。規格化は自動的に満たされている。
∫

1=
V

···

∫

PN (r0 , r00 , · · · , r(n) )dr0 dr00 · · · dr(n) .
(n)

V

17

(2.25)

第 2. 溶液状態を記述するための道具

第四分科会・佐藤啓文

しかしながら、困ったことにこの分布関数は、熱力学極限（密度 ρ ≡ N/V を一定にして、粒子数 N と体積 V について、N → ∞ と V → ∞ の極限）で、ゼロになってしまう。このことは、非常に大きな空間の中の微小体積要素に、非常に多くの粒子のうち特定のものがある確率を考えているのだから、殆ど自明である。熱力学極限操作は系のアンサンブル依(n)
存性を取り除く上で極めて重要であり、その意味ではこのままの PN を液体論の定式化に用いることはできない。そこで別の分布関数（generic distribution function）を導入する。
n 個の微小体積要素（dr0 , dr00 , · · · , dr(n) ）を考えるところまでは同じだが、それぞれの体積要素に入る粒子はどれでもいいとして確率を考える。すなわち、微小体積要素 dr0 に一個の粒子を入れる仕方は、どの粒子を選んでもよいので、全部で N 通りである。次に dr0 以
外の体積要素、例えば dr00 に残っている粒子から一つを選んで入れる仕方は (N − 1) 通りである。この作業を n 個の微小体積要素が埋まるまで続ければ、結局その仕方はN (N − 1) · · · (N − n + 1) = N !/(N − n)!,(2.26)

通りである。新しい分布関数の定義は

ρN (r0 , r00 , · · · , r(n) ) =
(n)

N!
(n)
P (r0 , r00 , · · · , r(n) ),(N − n)! N(2.27)である。この関数は熱力学極限で有限の値を持つ。最も簡単な相互作用がない場合（U = 0）を考えてみると、lim ρN (r0 , r00 , · · · , r(n) ) =(n)

V,N →∞

=

N!

lim

1

V,N →∞ (N − n)! V n

lim ρnV,N →∞N!
= ρn ,n
N (N − n)!

(2.28)(1)となる。一様な流体中では一体分布関数 ρN (r) は位置に依存しない定数（平均数密度）となることがわかる。
(1)

ρN (r) =N= ρ (constant).
V(2.29)

再びデルタ関数による表式を考えてみる。ある瞬間に空間中の微小体積素 dr に N 個のうちのどれでもよい、少なくとも一個の粒子を見いだす確率はN∑(1)
δ(r − ri )dr,νN (r)dr ≡i=1(2.30)

である。同様に n 個の微小体積要素の組 dr0 , dr00 , · · · dr(n) 中に N 個のうちのどれでもいい
n 個の粒子を見いだす確率は

νN (r0 , r00 , · · · , r(n) )dr0 dr00 · · · dr(n)
(n)

≡

N ∑
N
∑
i=1 j6=i···N∑δ(r0 − qi )δ(r00 − qj ) · · · δ(r(n) − qk )dr0 dr00 · · · dr(n) ,k6=i,j,···18(2.31)

第 2. 溶液状態を記述するための道具

第四分科会・佐藤啓文となる。この確率密度の熱平均値は既に導いた generic distribution function を与えることが確かめられる。
∫(1)< νN (r) >

···

=
V

∫

= N
V

∫ ∑
N
V i=1

···

∫

V

δ(r − ri )PN (r1 , r2 , · · · rN )dr1 r2 · · · drN

PN (r, r2 , · · · rN )dr2 r3 · · · drN(1)= ρN (r).

(2.32)

二体分布関数も同様にして得られる。

< νN (r, r0 ) >
(2)

∫ ∫

=
V

V

···

∫ ∑∑
V

= N (N − 1)

i

∫ ∫

j6=i

V

V

δ(r − ri )δ(r0 − rj )PN (r1 , r2 , r3 , · · · , rN )dr1 dr2 · · · drN

···

∫
V

PN (r, r0 , r3 , · · · rN )dr3 dr2 · · · drN

∫
∫ ∫
N!
=
· · · PN (r, r0 , r3 , · · · rN )dr3 dr2 · · · drN(N − 2)! V V
V= ρN (r, r0 ).
(2)(2.33)

この関数で r と r0 が空間的に十分離れているときは、二つの粒子がお互いの影響を受け
(1)
(1)
(2)
なくなるので、ρN (r, r0 ) = ρN (r)ρN (r0 ) となることは容易に理解できる。逆にある程度小さな r ≡ |r − r0 | では二つの位置に相関があるので、その相関 g (2) を次のように定義する。

g (2) (r, r0 ) = ρN (r, r0 )/ρN (r)ρN (r0 )
(1)

(2)

(1)

(2.34)

これは「二体相関関数」と呼ばれ、液体論では中心的な役割を果たしている関数である。
一様で等方的な流体では、二体分布関数・相関関数はともに相対的な距離のみに依存する。

ρN (r, r0 ) = ρN (|r − r0 |) = ρN (r),
(2)

(2)

(2)

g (2) (|r − r0 |) = ρN (|r − r0 |)/ρ2 .
(2)(2.35)(2.36)

このように g (2) は一つの点からの動径方向の距離のみに依存する関数であるので、動径分布関数とも呼ばれ、g(r) と書かれる。これ以上詳細には触れないが、これは X 線や中性子
線回折などから観測することが出来る量であり、液体の平均的構造を反映する「構造関数」の意味を持っている。
この節を終わるにあたり、この二体相関関数と熱力学量との関係について述べる。一成
分 N 粒子系の内部エネルギーは、前に導入した 12-6-1 型の様な対ポテンシャル u(r) の近似の範囲では< E >= N <N∑p2< u(|ri − rj |) >,
>+
2m
i>j

19

(2.37)

第 2. 溶液状態を記述するための道具

第四分科会・佐藤啓文と書くことができる。右辺の第一項は (3/2)N/β に等しい。第二項は可能な組み合わせ（対）N (N − 1)/2 に対してN∑< u(|ri − rj |) >

i>j

=

1
N (N − 1) < u(|r1 − r2 |) >2∫NN
1 N (N − 1) dr u(|r1 − r2 |) exp[−βU (r )]∫=2
drN exp[−βU (rN )]∫∫1 N (N − 1) dr1 dr2 u(|r1 − r2 |)∫=N2
drN e−βU (r )∫drN −2 e−βU (r )
N

∫
1∫
=
dr1 dr2 ρ(2) (r1 , r2 )u(|r1 − r2 |).
2(2.38)

さらに積分変数を (r1 , r2 ) から (r1 − r2 , r2 ) に変数変換すると、∫

V ρ2drg(|r1 − r2 |)u(|r1 − r2 |)
2 ∫
∫
1
1
=
N ρ drg(r)u(r) = N ρ 4πr2 drg(r)u(r).
22(2.39)

この式は直感的に非常に理解しやすい形をしている。すなわち、ある粒子から見たときに半径 r で幅が dr の球殻の中には g(r)dr = 4πr2 g(r)dr 個の粒子があり、その粒子との相互作用エネルギーは u(r) である。ポテンシャルエネルギー全体は全空間（r）に渡って積分することで得られる。系の内部エネルギーは
∫
<E>
31 1
=
+ ρ drg(r)u(r),N
2β 2(2.40)となる。この結果から想像できるように、系の熱力学的性質を個々の粒子と結びつけて理解するためには、動径分布関数を求める理論が必要である。

2.2.2

密度のゆらぎと OZ 方程式

この節では、動径分布関数を求めるための基本方程式である、Ornstein - Zernike 方程式を導く。
密度の空間的ゆらぎについて考える。密度場の平均値からのず̇れ̇は次式で定義される。

δρ(r) ≡ ν (1) (r) − ρ =N∑δ(r − ri ) − ρ.

(2.41)

i=1

一様な流体では一体分布関数は位置によらずに平均密度（ρ = N/V ）に等しいので、密度ゆらぎの一次のモーメントは常にゼロになる。

< δρ(r) >=< ν (1) (r) > −ρ = 0.

(2.42)
20

第 2. 溶液状態を記述するための道具

第四分科会・佐藤啓文

したがって液体中の密度の空間的ゆらぎを特徴づけるためには、二次以上のモーメントが必要になってくる。液体論において重要な役割を果たすのは二次のモーメント、すなわち
密度ゆらぎの相関関数である。
0< δρ(r)δρ(r ) > =*N∑δ(r − ri )

N
∑

+

0

δ(r − rj ) − 2ρ

j=1
2

i=1
0

*N
∑

+

δ(r − ri ) + ρ2

i=10

= ρδ(r − r ) + ρ g(r, r ) − ρ

2= ρδ(r − r0 ) + ρ2 h(r, r0 )(2.43)

最後の式で、第一項は i = j から来た寄与であり、粒子の自己相関に対応している。また、最後に定義した関数 h = g − 1 は、全く相関がなくなったときにゼロになる関数であり、全相関関数 (total correlation function) と呼ばれる。
前にも述べたように二体相関関数は X 線散乱や中性子回折によって直接観測可能な物理量であり、液体構造をそのゆらぎで表現する関数である。さらに、液体の微視的構造と熱力学量を結びつける最も基本的な関数でもある。二体相関関数については Ornstein-Zernike(OZ)
方程式と呼ばれる次のような関係式が知られている。
0∫0h(r, r ) = c(r, r ) +dr00 c(r, r00 )ρ(r00 )h(r00 , r0 ).

(2.44)

c は直接相関関数 (direct correlation function) と呼ばれる量である。OZ 式は Ornstein とZernike が半ば直感的に書き下した式であるが、その後多くの研究者の努力によって分配関数から導かれることが分かり、その性質も明らかになっている。
OZ 式の導出に先だって、グランドカノニカル統計集団（grand canonical ensemble; 大正準集合）を導入しておく。これまで、分布関数や相関関数を定義するために、NVT 一定の統計集団（正準集合）を扱ってきたが、この統計集団は粒子数のゆらぎが本質的である
ような問題には適応できない。そこで、ゆらぎを許すような統計集団を考える必要があるからである。グランドカノニカルの統計集団の定式化は次式で定義される分配関数が基礎
になる。
∫
∞
∑
zN ∫
Ξ =
· · · exp (−βUN ) dr1 · · · drNV
N =0 N ! V
=

∞
∑
zN
N =0 N !

ZN .

(2.45)

ここで z は活動度（activity）であり、

z = Λ−3 exp(βµ).

(2.46)と表される。グランドカノニカル統計集団の分布関数は以下の様に定義される。

ρ(n) (r0 , r00 , · · · , r(n) ) =< νN (r0 , r00 , · · · , r(n) ) >
∫
∫
∞
zN
1 ∑
N
· · · drN e−βU (r )=
Ξ N =0 N !
(n)

×
=

∑∑

···

i j6=i∞N∑∑δ(r0 − ri )δ(r00 − rj ) · · · δ(r(n) − rk )k6=i,j...

1z(n)
ZN ρN (r0 , · · · , r(n) ),Ξ N =0 N !

(2.47)
21

第 2. 溶液状態を記述するための道具

第四分科会・佐藤啓文(n)すなわち、カノニカル統計集団で定義した ρN を粒子数 N に関して平均したものになる。
ρ(n) の規格化は、
∫
V

···

*

∫

ρ
V

(n)

0

(r , · · · , r(n)0)dr · · · dr

(n)

=

+

N!
.
(N − n)!

(2.48)である。
さて、グランドカノニカルの分配関数から出発して OZ 式を導こう。UN は、i 番目の粒子と外場との相互作用 ψ(ri ) を含む一般的な形式で定義しておく。
N∑

UN (r1 · · · rN ) =ψ(ri ) + U 0 N (r1 · · · rN ).

(2.49)

i=1

一般化した活動度を z(r) = z exp{−βψ(r)} と定義すると、分配関数は

Ξ =

∫
∞
∑
1
N =0 N !

V

···

∫ ∏
N
V

z(ri ) exp(−βU 0N )dr1 · · · drN ,(2.50)iと書き直すことができる。これを z(r) で汎関数微分すると（Appendix 参照）、

δΞ[z]
=
δz(r)
=
=

∑ 1 ∫
N =0 N !

V

···

N ∫
∑ 1 ∑
N =0 N ! i=1 V∫ ∑δz(ri ) ∏V i=1 δz(r) j6=i

···

∫
V

δ(r − ri )

N ∫
∑

∑

z(ri ) exp(−βU 0N )dr1 · · · drN∏z(ri ) exp(−βU 0N )dr1 · · · drN

j6=i

∫

∏
1
···
z(ri ) exp(−βU 0N )dr2 · · · drN ,
(N
−
1)!
V
V
i−1
j=2
N =0

(2.51)となる。ここから以下の式が得られる。

ρ(r) =

δ ln Ξ
,δ ln z(r)ρ(2) (r, r0 ) =(2.52)

δρ(r).
δ ln z(r0 )(2.53)式 (2.53) の逆関係から直接相関関数 c(r, r0 ) を次式のように 定義 する。

δ ln z(r)δ(r − r0 )=− c(r, r0 ).
δρ(r0 )
δ(r)(2.54)(2.53) と (2.54) を汎関数微分の鎖則 (chain rule)、∫δρ(r) δ ln z(r00 )= δ(r − r0 ),δ ln z(r00 ) δρ(r0 )(2.55)に代入することで、OZ 式が得られる。OZ 式の本質は「全体」を表す h が「直接」の相互作用を表す関数 c の畳み込み積分の形で表現されているところにある。

22

第 2. 溶液状態を記述するための道具

第四分科会・佐藤啓文

この式は式 (2.54) で分かるように、いわば直接相関関数の定義を与えている式であり、二つの未知関数 h、c を含んでいるのでこれを解くためにはもう一つの関係式が必要である（積分方程式を完結するので、“closure relation”と呼ばれている）。代表的な関係式 として
HNC 近似がある。

c(r, r0 ) = exp[−βu(r, r0 ) + h(r, r0 ) − c(r, r0 )] − 1 − {h(r, r0 ) − c(r, r0 )}.

(2.56)

実際の手続きとしては、OZ 式とこれら closure relation を繰り返し計算によって解けば関
数 h、c を得ることができる。Figure 2.3 に液体アルゴンの動径分布関数を示した。HNC 近

図 2.3: 液体アルゴンの動径分布関数「液体の構造と性質」（戸田・松田・樋渡・和達、岩波書店刊）より抜粋
もとのデータは A.A.Khan, J. Chem. Phys. 134, A367, (1964)

似による計算は実験結果を非常によく再現していることが分かる。この OZ 式を分子（形を持った粒子）に一般化したものが RISM である。
HNC 近似の導出やそれ以外の closure relation、RISM の導出は分科会としての消化できる量を明らかに越えるので、ここではその理論の詳細には立ち入らないことにして結果を眺めることにとどめておく（[4]）。

23第 3 章 連続体モデル
連続誘電体モデルは、極性溶媒中では静電的な相互作用が支配的であることから、溶媒を分子として扱う代わりに分極可能の連続媒体と考える近似である。古くは Onsager やKirkwood のモデルまで遡る、いわば伝統的手法であり、量子化学計算との結合も 1970 年
代から進められていた。ここでは最近までに確立した幾つかの代表的方法を取り上げて、その方法論的側面について述べる。
連続誘電体モデルで扱われている「溶媒」は、シミュレーション等で用いられている分子構造を持ったモデルと比べると明らかなように、個々の溶媒分子を考慮しない、思い切った近似である。しかしながら溶液内分子の電子状態についての取り扱いは、この枠組みの中でさえ未解決の問題が残っており、代表的なモデルとしての意義は非常に深い。なお、
連続誘電体モデルはより一般的には連続体モデルの一部と見做すべきである。実際の問題
では、溶媒効果は極性のより強い溶媒で重要になると予想される。つまり静電的な効果が他の効果に比べて非常に大きいと考えるわけである。非極性溶媒はそれ以外の効果（分散力など）が支配的になると考えられるので連続体の、誘電体としての性質のみを考慮するだけでは不十分となってくる。しかし、ここでは連続体モデルの一般論を展開するのではなく、誘電体モデルについてだけを考えることにする。また、量子化学的ではない取り扱
いで、蛋白質などの大型分子を中心に溶媒和エネルギーを Poisson-Laplace 方程式を直接解くことで数値的に見積もる研究例は数多く知られているが、それらは他の総説 [2] 等にゆずることにする。

3.1

連続誘電体の一般的記述

ここでは注目する溶質分子が適当な大きさ・形を持つ連続誘電体中の空孔 (cavity) の中に存在していると考えよう。この時、溶質分子の電荷分布 ρsolute に基づいて周辺の溶媒（誘電体）が分極し、この変化がフィードバックして溶質分子の電子構造を変化させる。さらにこの電子構造変化は誘電体の分極に再び影響するので、全体として双方についての自己
無撞着的な取り扱いが必要となることが分かる。
連続体を考慮した全体の電荷分布を決定するためには、古典電磁気学に基づき、Poisson- Laplace 方程式を考える [3]。溶質分子上の電荷分布 ρsolute は原子核からの寄与 ρnuc と電子からの寄与 ρelec の和で表すことが出来、空孔内部の静電場 Φin と外部の静電場 Φout は以下の方程式を解けば得られる。

∇2 Φin (r)∇2 Φout (r)= −4πρsolute (r) = −4π[ρnuc (r) + ρelec (r)],
= 0.
24

(3.1)

第 3. 連続体モデル

第四分科会・佐藤啓文

Φout

Φout

AAA
AAA
Φin

AAAA
AAAAA
AAAA
AAAA
AAAAA
AAAA
AA
AAAA
AAAAA
AAAA
AA
AAAAA
AAAA
AA
AAAAA
AAAAA
ΦinAAAAA
AAAAA
AAAAA
AAAAA
AAAAA

AA
AA
AA

図 3.1: 連続誘電体モデル
分子形状の空孔 (PCM など) と楕円球の空孔 (SCRF など)

この際、境界条件として

Φin (r∗ )∂Φin (r∗ )∂n= Φout (r∗ ),∂= ε Φout (r∗ ).
∂n

(3.2)

が採用される。ここで n は空孔の表面に対する法線ベクトルを、r∗ は空孔と誘電体の境界面を、ε は連続体の誘電率をそれぞれ表している。空孔内部の静電場 Φin はさらに、直接溶
質分子の電荷分布がつくる部分 Φdirect と周辺の誘電体からの寄与による部分 Φreact に分割することができる。

Φin = Φdirect + Φreact .

(3.3)

右辺第二項の Φreact は反作用場と呼ばれる。溶質分子と周辺の誘電体との静電的な相互作
用エネルギー W は、反作用場と電荷分布 ρsolute との積を空孔内の体積について積分する
ことで得られる。
∫

W =r∈cavity 内部

drρsolute (r)Φreact (r).

(3.4)

従って溶媒和による自由エネルギー変化は

1∆Gsol = W,
2

(3.5)
25

第 3. 連続体モデル

第四分科会・佐藤啓文となる。
これを量子化学の手法と組み合わせるためには、反作用場が溶質分子の電子状態にも依存することに注意して、誘電体を含んだ新たな系のエネルギー Etot を定義すればよい。これは、孤立分子の電子ハミルトニアンを H、波動関数を Ψ と置くと[H + Hint ]Ψ = Etot Ψ,(3.6)と表すことができる。このうち Hint は先の相互作用エネルギー W と関連づけられる演算子である。

< Ψ|Hint |Ψ >= WM S .

(3.7)

従ってこのエネルギーに対する Fock 行列要素は(0)Fµν = Fµν
+ < µ|Hint |ν >,(3.8)となる。ここで F (0) は孤立分子に対する Fock 演算子である。
実際の問題を解く際には、新しく導入された誘電体による「溶媒和」を記述する演算子の表式を決定しなければならない。これを決定する二つの重要な要素は溶質分子の電荷分布の表現と空孔の形状・大きさである。現在までに提案されている主な方法は、おおよそ
Table 3.1 の様にまとめることができる。
表 3.1: 連続誘電体近似に基礎をおく手法

空孔の
形状
SCRF:
MPE:
COSMO:

球・楕円球
分子の形状

多極子展開
SCRF
MPE

self consistent reaction ﬁeld
multipolar expansion
conductor-like screening model

電荷分布の表現
相互作用点
連続分布
IC
—
GB/COSMO PCM
IC:
GB:
PCM:

image charge approximation
generalized Born methodpolarizable continuum model

理想的空孔の形状は、溶媒分子が入って来ることができない空間に相当していて、且つ溶質分子の全電荷分布を内部に含むようなものであろう。分子の形状にあわせた空孔は、この点のみに限定すれば最も合理的な選択の一つといえる。しかし、一方で空孔の形状を球または楕円球にすると、相互作用演算子 Hint が解析的に簡単に書き下すことが出来、溶質分子の核座標に対するエネルギー勾配法等の開発にあたって非常に有利になる側面もある。
溶質分子の電荷分布の表現に関しての二つの方法は既に 2.1 節で説明した通りである。多重極子展開は伝統的かつ典型的な方法の一つであるが、一方で、代表的なシミュレーショ
ン用パラメーターセットの殆どは相互作用点モデルの考え方に基づいている。すなわち、分子上で幾つかの相互作用点を選び、これらの上での電荷分布（殆どの場合は点電荷）の和で全体の電荷分布を表現する方法である（式 2.6）。とりわけ、分子の形状が複雑になると多重極子展開は低次まででは収斂しにくくなることが知られており、相互作用点を用いる方が有利であると考えられる。

26

第 3. 連続体モデル

第四分科会・佐藤啓文

3.2

連続誘電体モデルの例

3.2.1

Self Consistent Reaction Field(SCRF) 法

この方法は最も単純で最も古い方法である。半径が a である球形の空孔を考えると、その中心に多重極子を置いたときの反作用場は、Legendre 多項式を使って、次の様に解析的表現で表すことができる。
∞∑

l(l + 1)(ε − 1) rl+1 ∑Φreact =−Blm Ylm (θ、φ).
2l+1
l + (l + 1)ε a
l=0
m=−l

(3.9)

このうち最も単純で重要な例は、l = 0 の項、すなわち空孔の中心にある一個の点電荷だけを考えたケースであり、その溶媒和自由エネルギーは Born の式として知られている。

∆G = −

ε − 1 q2,2ε a(3.10)

ここで、q は空孔内にある電荷を示す。電気的に中性な化学種ではこの項は消えてしまい、
有意の最低次は l = 1 つまり双極子モーメント (µ) からの寄与になる。

∆G = −

ε − 1 µ2.
2ε + 1 a3(3.11)

量子化学的計算への最初の応用は、Rivail 達によって 1973 年に報告された、CNDO の例であろう。彼らは l = 7 までの多重極子を考慮して、その収束性や自由エネルギーについて議論した [28]。彼らの Fock 演算子は簡単には(0)Fµν = Fµν+

∑ ∑l、l0 m、m000、m
Mlm flm< µ|Mlm
0 |ν >,、l0(3.12)0、mと書ける。ここで Mlm は多重極子演算子であり、flm
は空孔の大きさと誘電率によって、l0決まる定数である。ab initio 法への応用は 1983 年に行われ、この後さらに楕円球型の空孔の場合や、1990 年代に入ってからは多体摂動論の考慮、エネルギー勾配法の開発や、多中
心多極子展開の方法への拡張を行っている [29]。同様のハミルトニアンを基礎として最初に ab initio 法へ応用した例は 1974 年に Hylton、 Christoﬀersen、Hall によって報告された[30]。
SCRF 法の枠組みでもう一つ特筆すべきは Wong、Frish そして Wiberg [31] による報告であろう。彼らは双極子の項までを入れた（中性分子の場合は双極子の項のみになる）Fock
演算子、

(0)

Fij = Fij −ε−1 µ< i|µ̂|j >,2ε + 1 a3(3.13)

を基に、二次までのエネルギー勾配法など様々な計算を行った。この方法は 最もポピュラーなプログラムパッケージの一つである Gaussian シリーズにリンクされ、電子相関を考慮した方法 [32] などの開発・拡張が行われ、広く利用されている。現在はかなり多くのプログ
ラムパッケージが、この方法を持っているようである。

27

第 3. 連続体モデル

第四分科会・佐藤啓文

SCRF 法の最大の利点は簡便さにあるといえる。そのため、エネルギー勾配法の開発やCI 法、MCSCF 法、あるいは密度汎関数理論といった他の方法への拡張も比較的容易に行われてきた。しかしながら、モデルの近似があまりにも粗野なので定量的な議論に耐えうるかどうかは疑問が残ると言わざるを得ない。数多くの指摘がなされている様に、一般に多重極子展開の収斂は遅く、低次での打ち切りがどの程度の誤差を与えるかを注意深く慎重に調べる必要があろう。また空孔の形状が注目する分子の形状とは大きくかけ離れてしまう上に、球形の場合ならば半径をどの様に決定するか、加えて楕円球の場合は軸の方向を如何に取るか等、任意性の残り易い幾つかの問題点があることを指摘しておく。

3.2.2

Generalized Born(GB) 法

Cramer と Truhlar によって 1989 年に提案された [33] この方法では、分子全体の溶媒和
エネルギーが、相互作用点（溶質分子の各原子）上の有効電荷に対する溶媒和エネルギーの和になると仮定している。各相互作用点に於ける溶媒和エネルギーを Born の式で見積もる。一見すると、多中心展開の SCRF 法で、各点上での展開を第一項目で切ったことと等価に思われるが、各相互作用点間の関係は経験的に扱われているので、この近似をさらに進めたことに相当すると考えられる。彼らは半経験的方法 (AM1) に基づいて、次のような Fock 行列を提案している。
(0)Fµν = Fµν
+ δµν

ε−1∑
qk γkl .
εk(3.14)

ここで、

qk = Zk −

∑

Pµν ,

(3.15)µ∈kであり、さらに、

γkk =

1
,
ak

γkl = √

1

2
2
2
rkl
+ ak al [exp(−rkl/dkl ak al + Ckl (rkl))],

(3.16)

である。式中の ak は相互作用点 k に於ける有効半径、rkl は相互作用点間の距離、dkl と Ckl
は誤差を補正するためのパラメーターと関数である。この方法では必要となる様々なパラメーターの決定法や、それらの物理的な背景がわかりにくいといった問題点が残っている
ものの、実験で観測されている溶媒和エネルギーの値を非常によく再現することが知られている。彼らはさらに非極性分子にも応用できように、分散力を取り扱うための提案、パラメータの改善などを精力的に行っており、SMx モデルとしてほぼ完成したパッケージになったようである。

28

第 3. 連続体モデル

3.2.3

第四分科会・佐藤啓文

Polarizable Continuum Model (PCM)

この方法は Tomasi らのグループにより 1981 年に発表され、Apparent Surface Charge(ASC) 法とも呼ばれていた。これまでに様々な応用・改良が行われており [34]、HartreeFock(HF) 法の他に UHF (Unrestricted HF)、MCSCF、MBPT、CASSCF 、MR-SDCI、
DFT 等についての報告がされている。最重要な特徴は分子形状の空孔を扱えることである。基本的方針として、反作用場を連続誘電体との境界面上に於ける（見かけ上の）電荷分布で表現することとして、そこから Hint を数値的方法で計算する。
境界面上の点 S に於ける電荷分布 σ(S) は次の式で計算できる。

σ(S) = −P · n =ε−1ε−1∇Φout · n、=∇Φin · n
4π
4πε

ε−1 ∂
(Φsolute + Φreact ),4πε ∂n∫σ(s) 2d s,
Φreact (r) =|r − s|=(3.17)

ここで P は分極ベクトルであり、数値微分でこれを求めている。Fock 行列要素は∫(0)Fµν = Fµν
+ < µ|

σ(s) 2d s|ν > .
|r − s|(3.18)と書き下すことが出来る。実際には、境界面を適当な大きさの面積素に分割して、積分をそれらの和で書き換える。

qk = ∆Sk σ(rk ).
∑qkΦreact (r) ∼.
k |r − sk |(3.19)

従って境界面を自由に設定することができ、通常は 各原子を中心とする van der Waals 半
径程度の大きさの球の集合体として空孔を表している。自由エネルギーの大きさは、他の方法と同じように空孔の形状や大きさに敏感であるが、球形の空孔を使っている場合に比べて実験値との比較に耐えうる範囲で再現できている様である。

3.2.4

Image Charge(IC) 法

導体の古典電磁気学で、適当な仮想電荷（鏡像電荷）を置くことで誘導された電場を代表させる方法がよく知られている。近似的にではあるが、この方法を誘電体の場合に拡張したのが Friedman であった [35]。例えば球状の空孔内に単一の電荷 (大きさ q) が空孔の中心から距離 r だけ離れた点にある時、反作用場による溶媒和自由エネルギーは次式で近似される。

∆G(0) (r) = −1
1 ε − 1 q2.
2 ε + 1 a 1 − r2 /a2(3.20)

前述の式 (3.10) と比較しても分かるように、これは r = 0 の点で Born の式に相当する（但し、(ε + 1)−1 のオーダーの誤差を含む）。この近似の大きな利点は、球内部の任意の位置

29

第 3. 連続体モデル

第四分科会・佐藤啓文にある電荷とその反作用場との相互作用が、空孔外の仮想電荷との相互作用で計算することができるところにあり、さらに空孔内部に電荷が複数存在する場合も、それら各々の仮想電荷を考えることで相互作用を計算できるところにある。球状の空孔に制限されるものの、多極子展開の全項を考慮できるので、注目する溶質分子の電荷分布を直接的に反映した反作用場を得ることができると言っても良いだろう。また、溶質分子の核座標に対する
エネルギー勾配法の表式が比較的簡単に書き下せるといった利点もある。
量子化学計算との結合は、最初 Karlström によって発表され [36]、後に筆者らによってエネルギー勾配法等を含めた研究がなされた [37]。

3.2.5

Conductor-like Screening Model(COSMO) 法

Klamt と Schüürmann は、PCM と IC を組み合わせた様なモデルを提案した [38] 。PCMと同じように空孔表面での有効電荷を求めることで任意の形状の空孔を扱える様に改善を行い、同時に境界条件として IC と同じように誘電率が ∞ のときに正しくなる近似を導入する。彼らの結果は古典的点電荷の取り扱いと半経験的手法との組み合わせだけであったが、のちに Troung らのグループによって非経験的手法への拡張がなされている。
PCM と同じ様に、境界表面を M 個の小さな面積要素（面積 Sµ ）に分割する。それぞれの電荷密度を σµ とすると、任意の点 ri に単位電荷が作るポテンシャル biµ は
∫

biµ =

∫
σµ
1
1
d2 r =d2 r|Sµ | Sµ |r − ri |Sµ |r − ri |≈ |r − ri |−1 .

(3.21)(3.22)

である。同じ様に表面電荷分布どうしの相互作用は∫

1d2 rd2 r0 ,0
Sµ Sν |r − r |(3.23)µ 6= ν のとき aµν ≈ |rν − rµ |−1 ,(3.24)µ = ν のとき aµµ ≈ 3.8 |Sµ |,(3.25)1aµν =|Sµ ||Sν |
この要素は

√

である。これらの行列要素と溶質分子上の有効電荷（Q）を用いて、相互作用エネルギー∆E を最小にする表面電荷の組を求めると、最終的に

1∆E = − QBA−1 BQ,
2(3.26)

を得る。面積分と相互作用の行列要素について大胆な近似を導入しているものの、PCM にかなり近い結果を与えているようである。

30

第 3. 連続体モデル

3.2.6

第四分科会・佐藤啓文

連続誘電体近似の応用例

最後に化学反応への適用例の一つとして、筆者が行った鏡像近似法と ab initio との結合と、その応用を示す [37]。分子系の溶媒和した自由エネルギー G は、溶質分子の持つ Esolute
ˆ でと溶質分子の電荷分布 (ρnuc (r)、 ρelec (r)) とそれぞれの鏡像電荷を生成する演算子 (Im示す）を用いて次の様にクーロン相互作用（El）の和で表せる。

1
ˆ nuc } + El{ρelec Imρˆ nuc }
[El{ρnuc Imρ
2ˆ elec } + El{ρelec Imρˆ elec }] + Vex .
+ El{ρnuc Imρ

G = Esolute +(3.27)

Vex は溶媒分子の電子との交換相互作用による排斥ポテンシャルである。Hartree-Fock の枠組みでは、これから Fock 行列要素が密度行列 D を使って次のように導かれる。
0Fµν = Hµν+∑

1(µλ|σν)]
2
∑
1
1
ˆ
ˆ
+ (λσ|Im|µν)].
+Dλσ [ (µν|Im|λσ)22λσDλσ [(µν|λσ) −λσ(3.28)

右辺第一項目は「原子核と電子の鏡像」、「電子と原子核の鏡像」、交換相互作用を含んだ
一電子積分を示し、末項が「電子とその鏡像」の相互作用を示す二電子積分を表している。
鏡像との相互作用に含まれている二体の演算子をさらに次のように書き直す。
(∞∑
1
r2 − r Q
11==|r1 − r2 ||(r1 − rQ ) − (r2 − rQ )||r1 − rQ | l=0 r1 − rQ)l

Pl (cosθ).

(3.29)Q は二つの原始関数の中心位置を示す。原始関数 Gi 、 Gj および Gk 、 Gl に対する二電子
積分の項は
ˆ k Gl >< Gi Gj |Im|G
=

∞
∑
l=0

Pl (cosθ) < Gi |1ll+1 |Gj >< Gk |(r2 − rQ̂ ) |Gl > .
|r1 − rQ̂ |(3.30)となり、この l についての展開の第一項目のみをとれば一電子積分で近似できるので、大幅な計算時間短縮が可能になる。一電子積分についても同様の近似が可能である。
SN 2 反応の結果を示そう。この反応はよく知られた有機反応であるが、同時に溶媒和の影響が極めて顕著な例としても知られており、これまでにも多くの量子化学的研究が発表されている。ここでは次の様な最も典型的反応例について取り上げる。

Cl− + CH3 Cl → ClCH3 + Cl−(3.31)−
反応座標 s として、二つの Cl− –C の距離の差 s = r(Cl−1 −C) − r(Cl2 −C) を取ると、対称的な s = 0 なる点が遷移状態となる。溶媒としては水（ε = 80.0）を仮定している。
この s を横軸にとった反応のポテンシャル面を Figure 3.2 に示す。横軸は反応の対称性から半分だけ（遷移状態→生成系、あるいは、遷移状態←反応原系）を、縦軸は気相中・溶液中それぞれの遷移状態でのエネルギー値を 0.0 kcal/mol とした相対値で示してある。

31

第 3. 連続体モデル

第四分科会・佐藤啓文

Relative Energy / kcal mol - 1

0
in gas phasein aqueous solution
in aqueous solution(with integral approx.)

-5

-10

-15

-20

0

1
2
3
4
5
reaction coordinate s / bohr6

図 3.2: SN 2 反応のポテンシャル面

実際は系全体で電荷を持っているので、溶液中では Born の式で見積もられるのとほぼ同
じ 24 kcal/mol の気相中に対しての安定化エネルギーを得ている。
気相中では、これまでよく知られているように、s=3.0[bohr] 付近でイオン–双極子の安定錯体が存在していることがわかる。一方、溶液中ではこれらの錯体が消失しており、反応の活性化エネルギーが気相中の値（この錯体からで 15. 5 kcal/mol）に比べて約 20 kcal/molと高くなっていることが分かる。また、先の分子積分に対しての二つの方法を図中に示してあるが、実在の電荷分布のその仮想電荷の距離が十分遠いため、近似が極めて精度の良いものになっていることが分かる。
これまでの実験・シミュレーションで得られている水中での活性化エネルギーの値は約26 kcal/mol であり [41]、誘電体近似で得た値は、これに近いものの、若干小さい。数値の差違の原因としては、球状の空孔を用いているために排除体積が過大に評価されている点、また解離極限での値を知ることができない点などが考えられる。しかしながら、連続誘電体モデルの計算は非常に簡便であるにも拘わらず、少なくとも定性的範囲では満足いく結果が得られたことも強調しておきたい。

32

第 3. 連続体モデル

3.3

第四分科会・佐藤啓文

連続誘電体近似のまとめ以上挙げた方法の他に代表的な方法として、Hoshi らによる matrix inversion 法 [40] などが挙げられる。これは反作用場の離散的表示に基礎を置く点で PCM の考え方に極めて近いが、より解析的な表記になっている。ここまで概観してきた幾つかの手法の中で、近似的方法としてどれが最も優れているかを決定するのは難しいが、理論的経緯も明確で溶媒和エネルギーの絶対値も求められている点で、PCM 法が有力な候補として挙げられるだろう。
すぐ気づくことだが、溶媒を連続誘電体としている以上、空孔の形状・大きさが決定す
れば、Poisson-Laplace 方程式の解は一意に決まる。この式を直接数値的方法で解ければ1 、それが連続誘電体近似という枠組みの中で得られる最良の答えであり、PCM はそれに最も近い値を与える近似方法である、ということである。
いずれの方法を用いた場合も、連続誘電体で近似していることから来る問題点を十分把握しておくべきである。例えば、空孔の大きさ（球の場合は半径）の決定には任意性が残り易く、しかもエネルギー等の大きさはこのパラメーターに対して敏感であることが知られている [34]。また溶媒の個々の分子構造を無視しているために水素結合に代表されるような局所的な相互作用の取り込みが難しい。
連続誘電体近似は量子化学計算との結合を通して定性的に分かり易い結果を与えて来た点で非常に意義深く、今後も第一段階の近似としての役割を果たしていくと考えられる。
しかしながら、より microscopic で定量的な情報を引き出していくためには困難が多いと言わざるを得ない。

1

実際、量子化学計算の直接的結合も試みられている。

33

第4章

QM/MM 法

これまでの議論からも分かるように、連続誘電体モデルの抱える問題の原因の一つは、溶媒の分子的描像とそれらが作り出す溶媒構造を極端に簡単化あるいは無視してしまったことにある。この直接的解決法は明らかである。全ての分子（溶媒も）を直接扱って、重要な部分を量子化学計算でそれ以外は古典論的手法で記述すればよい。溶媒の配置について
はシミュレーションを行う必要があり、莫大な時間が必要となるので、量子化学計算の部
分は多少粗い近似のもとで行う方が現実的であろう。この溶質の電子状態の記述について
は、これまで様々な提案がなされた。例えば、幾つかの電子状態間の遷移を取り扱う溶液
内電子移動反応の場合は、着目する電子励起状態で全体の波動関数を展開すればよい [42]。
また、結合の解離を記述するためには予め Valence Bond の描像を入れた半経験的な電子状態の記述をしておき、結合／解離の二つの状態で展開してやればよい（EVB 法）[43]。
QM/MM 法は名前から想像できるように量子化学（QM）計算と分子力場（MM）計算を混合した方法である。発想そのものが明解であり、1970 年代ころの比較的早い段階か
ら、孤立分子の π-σ を別々に取り扱ったり、生体分子でいわゆる活性中心とそれ以外を分離して扱うなど、様々な系に応用されてきた。しかし、溶液を扱った最初の例は、80 年代の Warshel や Karplus らによる研究のようである。一般にシミュレーションが非常に計算時間を要するので、計算機環境の整備が一つの大きなきっかけとなって、手法として現実的になってきたのであろう。生体分子内反応については QM と MM の接合部分の取り扱いが簡単ではないのだが、溶媒分子が反応に直接関与せずに場を提供しているタイプの溶液
内反応の場合はそのような困難はない。最近では J. Gao が精力的にこの方法に基づいた溶液内有機化学反応の研究を行っており、連続誘電体を越えていくための方法としてその重要性を増しつつある。
QM/MM 法では系のハミルトニアンを以下のように定義する。
0el
vdW
Hef f = HQM+ HQM/M
M + HQM/M M + HM M

(4.1)

0
は量子化学的に記述された（電子）ハミルトニアン、HM M は分子力場のハ
ここで、HQM
el
vdW
ミルトニアン、HQM/M
M と HQM/M M はそれらの結合部分に相当し、前者は静電的な相互作用、後者は van der Waals 相互作用を表している。実際には
el
HQM/M
M = −

vdW
HQM/M
M =

S ∑
N
∑
eqs

+

S ∑
M
∑
Zm qs

,

s=1 i=1 ris
s=1 m=1 Rms
[
(
(
)
) ]
S ∑
M
∑
σms 12
σms 6
4εms
−
.
Rms
Rms
s=1 m=1

(4.2)
(4.3)

である。ここで、N は溶質分子の電子の数、M と S は溶質、溶媒分子の数を表している。
式中のパラメーターは通常のシミュレーションで使うものを用いる。

34

第 4. QM/MM 法第四分科会・佐藤啓文

計算は、(1) 溶媒分子の配置を発生する。(2) 上記のハミルトニアンを使って分子軌道計算を行う。この手続きを繰り返せばよい。モンテカルロ法・分子動力学法いずれとも組み
合わせることができる。しかし、実際の計算は非常に時間がかかるので非経験的分子軌道計算による例は極めて少なく、大部分の研究は半経験的分子軌道法 (AM1 や PM3 による）
によるようである。
実際の計算結果については文献 [6] を参照することにする。

35

第5章

RISM-SCF/RISM-MCSCF 法

誘電体近似を超えて実際に計算可能である、もう一つの方法が、この章で取り上げる
RISM-SCF 法である。QM/MM 法は半経験的分子軌道法と組み合わせている事例が殆どであり、非経験的な方法との組み合わせという点では、現時点では RISM-SCF/MCSCF 法が、事実上唯一の方法である。
1993 年に Ten-no らによって提案された [45] この方法の特徴は、溶媒構造を RISM 方程式の援用で求めるところにある。後にこの方程式は筆者らによって再導出された。[46] こ
こでは後者の手続きに従って、この方程式を導く。

RISM-MCSCF 方程式

5.1

溶質分子に対しての溶媒分子の分布関数は、それらの分子間二体相関関数 (h )、直接相関関数 (c) と分子内相関関数 (ω) に関する積分方程式

ρhρ = ω ∗ c ∗ ω + ω ∗ c ∗ ρhρ,

(5.1)

を解くことで得られる [10]-[11]。ここで「*」は畳み込み積分を、ρ は系の平均密度を示し、各相関関数は原子対を添字とする行列である。実際に解くためには”closure relation”と呼ばれるもう一つの式が必要であり、代表的なものとして HNC(Hyper Netted Chain) 近似

cαs (r) = e[−βuαs (r)+hαs (r)−cαs (r)] − [hαs (r) − cαs (r)] − 1,

(5.2)などが挙げられる。式中 β = 1/kB T であり、kB 、T は Boltzmann 定数と温度を、また、α、s 等の添字は各々溶質・溶媒分子上の相互作用点を表している。
全系の Helmholtz 自由エネルギー A を、溶質分子起因のエネルギー Esolute と溶質・溶媒
分子間の相互作用による excess chemical potential ∆µ の和で定義する。

A = Esolute + ∆µ.

(5.3)

ここで

Esolute =∑ Za Zb
a>b

=

Rab

∑ Za Zb
a>b

Rab

+

∑

γij hij +1∑
Γijkl (φi φj |φk φl )2 ijklγij hij +1∑
Γijkl < φi |gkl |φj >,2 ijkl

ij

+

∑
ij

36

(5.4)

第 5. RISM-SCF/RISM-MCSCF 法

第四分科会・佐藤啓文であり、 γij と Γijkl は分子軌道と CI 係数から決まる vector coupling 係数である。他の記号は通常使われている意味である。一方、∆µ は、Singer と Chandler によって導かれた
(

)

ρ ∑∫
1
1
∆µ = −
dr cαs − h2αs + hαs cαs ,
β α,s
2
2

(5.5)

を採用する。この自由エネルギー A は溶媒構造を記述する相関関数 hαs 、cαs と tαs 、溶質の電子状態を表す MO 係数 vri と CI 係数 CI の汎関数と見なすことができる。そこで直交条件を課して Lagrangian L を次の様に定義する。

L = A[c, h, t, v, C] − E(∑)

CI2 − 1

−

∑∑
i

I

εim (Sim − δim ).

(5.6)

m

ここから変分的に∫{}ρ∑
dr e[−βuαs (r)+tαs (r)] − 1 − hαs (r) δtαs (r)δL = −β α,s+ {−tαs (r) + hαs (r) − cαs (r)} δhαs (r)+ −hαs (r) + ρ−1


+2

{
∑ ∑
I

+2

∑

∑




ωαγ ∗ cγs ∗ χss0 (r) δcαs (r)

γ,s0



}

HIJ CJ − ECI δCIJ< δφi |γij h +i,j∑Γijkl gkl

kl

(

)

ρ ∑∫
∂
− γij
−
dre[−βuαs (r)+tαs (r)] − εij |φj >,
bλ
∂qλ
β α,s
λ
∑

(5.7)

を得る。Singer と Chandler の論文に示されているように、この最初の三つの括弧の中は、それぞれ HNC closure、関数 tαs (r) の定義式、RISM 方程式になっている。次の項は通常の永年方程式で、最後の項が RISM-MCSCF 方程式である。この項に含まれている qλ は、
溶質分子上の相互作用点に於ける有効電荷であり、適当な密度演算子 bα を用いて波動関数と関係づけられる。

qα (v, C) =∑γij < φi |bα |φj > .

(5.8)

ij

さらに溶質分子上の点 λ に働く静電ポテンシャル (V)λ を要素とするベクトル V を導入す
れば、
(

)

∑ ∫ gλs (r)∂ρ ∑ ∫ [−βuαs (r)+tαs (r)]−eqsdr ≡ (V)λ ,dr = ρ
∂qλ
β α,s
r
s

(5.9)

であるから、Lagrange の未定乗数 εim を使って RISM-MCSCF 方程式の Fock 行列は、

εim =

∑
j

=

∑

< φi |γmj h +∑Γmjkl gkl − γmj Vt · b|φj >kl< φi |Fmj |φj >,(5.10)

j37

第 5. RISM-SCF/RISM-MCSCF 法

第四分科会・佐藤啓文

CI 行列要素は、
HIJ =∑γijIJ < φi |h − Vt · b|φj > +ij

1 ∑ IJΓ < φi |gkl |φj >,2 ijkl ijkl

(5.11)と決定される。HF 法の場合の Fock 行列要素は、Ten-no らによって示されたそれと一致する。
式の導出からも明らかなように、この方法では溶質の電子状態は溶媒分子の分布（動径分布関数）によって決まり、また動径分布関数は溶質分子の電子状態によって決定される
ので、双方が自己無撞着的に解かれる。このため、溶質の電子状態と溶媒分子の分布が、互いを反映した解を得られるのが特徴である。また、この枠組みで溶質分子の電子励起状態に対しての解も得られることが分かる。
また、ここでは省略するが、溶質分子の核座標に対してのエネルギー勾配法も、以上の枠組みの中で容易に導出することが出来る。

溶液内化学系への応用

5.2

この方法の応用例をいくつか紹介する。

5.2.1

メチルアミン類の塩基性順序

はじめに Kawata らによる、メチルアミン類の塩基性の順序の不規則性に関する研究を紹介する。水溶液中でのアンモニアあるいはメチルアミン類とその共役酸の解離反応の自由エネルギーを並べると、ジメチルアミン（(CH3 )2 NH）とトリメチルアミン（(CH3 )3 N）の間で、置換したメチル基の数から期待される順序との逆転が観測される。

NH3 > (CH3 )NH2 > (CH3 )2 NH < (CH3 )3 N

(5.12)

RISM-SCF 法で計算した自由エネルギー変化を Figure 5.1 に示す。Fig. ではアンモニ
アとその共役酸とでの解離エネルギーを基準とした相対値で示してある。ここから、順序の逆転は溶質分子からの寄与と溶質-溶媒間の相互作用からの寄与の微妙なバランス関係に起因していることが分かる。
また、ここには示さないが、メチルアミン類によって窒素と酸素間の動径分布関数は大きく変化しており、溶液内化学種を取り扱っていくためには、溶質分子の電子状態を正しく記述するだけでなく、溶媒構造をも適切に記述する事が重要であることが、この研究か
ら分かる。

5.2.2

ホルムアミドの互変異性化

ホルムアミドの互変異性化は分子内プロトン移動の最も小さな系としてとくに理論計算に於ては古くから研究対象になっている。

38

第 5. RISM-SCF/RISM-MCSCF 法

1

第四分科会・佐藤啓文

5

0

1

5

l

o

0

m

/

l

a

k

c

5

0

1

1



G

(

s

o

l

u

t



G

(

s

o

l

v

e



G

)

n

t

)

5



2

e

































0

0

0

.

.

5

1

2

1

N

u

m

b

e

5

r

o

2

.

5

f

C

H

3

3

図 5.1: メチルアミン類の自由エネルギー変化▲（溶質・溶媒間）、■（溶質）、●（全エネルギー）、○（実験値）
J. Am.Chem.Soc.、 117、 1638、Fig.2

NH2 -CHO *) NH=CH(OH)
この両辺の自由エネルギー差は溶媒によってどの程度変化するものであろうか。溶媒の個性は連続誘電体近似では ε の数値としてしか反映されない。また QM/MM 法では水以外の溶媒についての研究はまだそれほど進んでいないようである。
石田らは上の反応を様々な溶媒中で調べた [47]。
まず、溶媒の個性について計算によって得られた結果と実験事実との対応を見てみる。
測定結果に含まれる溶媒和効果を系統的に整理・解釈するための半経験的なパラメーター
として「Taft のパラメーター」がある。RISM-SCF/MCSCF 計算で得られた水素結合部分に相当する対相関関数（動径分布関数）のピークの高さの対数を取ったものと Taft のパラ
メーターの水素結合能力を示す指数とを比較したのが Figure 5.2(a)(b) である。この結果
は RISM-SCF/MCSCF 法が微視的な溶媒和構造をよく再現していることを示しており、同時にこれまで知られていた半経験的なパラメータに ab initio（第一原理からの）理論的手法でもって物理的意味づけをすることができた、とも言えるだろう。
Figure 5.2(c) に Onsager - Kirkwood のパラメータ (ε−1)/(2ε+1) と RISM-SCF/MCSCF法で得られた溶媒和エネルギーをプロットしてある。図での線型関係からのずれが、いわば、誘電体近似では表せない溶媒の個性を表していると解釈できる。とくにアセトニトリ
ルのずれは非常に大きいことが分かる。

39

第 5. RISM-SCF/RISM-MCSCF 法

第四分科会・佐藤啓文

図 5.2: 様々な溶媒中でのホルムアミドの互変異性化反応40

第 6 章 まとめと展望
本稿では溶液内分子の電子状態を取り扱うための手法について概観してきた。繰り返しになるが、多原子系である分子の問題を扱う上で、鍵となるのは溶質・溶媒分子を含めた
全系の持っている莫大な数の自由度と、それぞれの分子（または分子間）の量子力学的な側面を如何に整理して扱って行くかという点であろう。最終的に前者を扱うためには統計力学的な方法を導入することが不可避であった。一方、後者は電子ハミルトニアンの近似といった形で、より効率の良い方法が模索され続けている。
溶質分子の電子構造と溶媒構造の自己無撞着性という観点からは、現段階でも、QM/MM
法と RISM-SCF/MCSCF 法が実際の問題に対して有用な情報を提供してくれることが期待される。しかしながら、これまでの議論からも明らかなように、次の新しい方法論による解決を待つ幾つかの問題点が残っている。例えば、溶質・溶媒分子間に働く力の中で最も支配的なクーロン的相互作用を問題として、他の相互作用はパラメーターに暗に含まれているとするものが多かった。しかし、より一般的には二つの分子間に働く、これ以外の交換、電子移動や溶媒側の分極といった相互作用を扱っていく必要があろう。勿論、分子間の相互作用の記述がそれほど簡単ではないことは本論で触れた通りである。また、 溶媒分子の寄与を平均場的に取り込む考え方は、扱う問題によっては明らかに破綻することも指摘しておく。
溶液内の電子状態を調べる本格的な試みは、まだ始まったばかりである。しかしながら、本稿で紹介した幾つかの手法の確立によって、溶質・溶媒の構造・反応に対しての様々な情報を引き出すことが可能となったわけである。さらなる発展によって、従来の気相中電子状態理論に匹敵する有効性を有す理論体系が構築されることを期待し、また目指したいと思っている。

謝辞
分子研・平田文男教授には、日頃の討論ばかりでなく、本稿作成にあっても貴重な意見を頂き、Appendix の転載を快諾して頂いた。また、高山敦子氏には原稿の査読をしていた
だいた。この場を借りてお礼申し上げます。

41

関連図書
[1] J. Tomasi and M. Persico, Chem. Rev., 94, 2027, (1994). （量子化学との結合を志向
した連続誘電体について）[2] M. E. Davis and A. McCammon, Chem. Rev. , 90, 509, (1990). （古典的取り扱いの連続誘電体について）[3] C. J. F. Böttcher, 1973, Theory of Electric Polarization (Elsevier).（誘電体を中心とする教科書）[4] P. J. Rossky, Ann. Rev. Phys. Chem., 36, 321, (1985). （RISM についての総説）[5] F. Hirata, Bull. Chem. Soc. Jpn., in press, (1998) （RISM についての総説）[6] J. Gao, Acc. Chem. Res., 29, 298, (1996). （QM/MM 法についての総説）[7] 佐藤啓文、物性研究, 68, 434, (1997).
[8] J. P. Hansen, and I. R. McDonald, Theory of Simple Liquids, Academic press, 1986.
[9] 平田文男、湘南レクチャー「溶液の構造と緩和」講義資料[10] H. C. Andersen and D. Chandler, J. Chem. Phys. 57, 1918, (1972) ; D. Chandler
and H. C. Andersen, J. Chem. Phys. 57, 1930, (1972).
[11] F. Hirata and P. J. Rossky, Chem. Phys. Lett. , 83, 329, (1981) ; F. Hirata, B. M.
Pettitt and P. J. Rossky, J. Chem. Phys. , 77, 509, (1982) ; F. Hirata, P. J. Rossky
and B. M. Pettitt, J. Chem. Phys. , 78, 4133, (1983).
[12] J. T. Blair, K. Krogh-Jespersen and R. M. Levy, J. Am. Chem. Soc. , 111, 6948,(1989).
[13] H. Fukunaga and K. Morokuma, J. Phys. Chem. , 97, 59. (1993).
[14] J. Zeng, J. S. Craw, H. S. Hush and J. R. Reimers, Chem. Phys. Lett. , 206, 323,
(1993) ; J. Zeng, J. S. Craw, H. S. Hush and J. R. Reimers, J. Chem. Phys. , 99,1482,1496, 1508,(1993).
[15] J. A. C. Rullmann and P. T. van Duijnen, Mol. Phys. , 61, 293, (1987).

42

第 6. まとめと展望

第四分科会・佐藤啓文[16] M. F. Ruiz-López, D. Rinaldi, and J. Bertrán, J. Chem. Phys., 103, 9249, (1995) ;R. Bianco, S. Miertuš, M. Persico, and J. Tomasi, Chem. Phys. , 168, 281, (1992).
[17] D. Chandler and P. G. Wolynes, J. Chem. Phys., 74, 4078, (1981) ; R. N. Barnett,
V. Landman, and A. Nitzan, J. Chem. Phys. , 89, 2242, (1988) ; J. Schnitker and P.
J. Rossky, J. Chem. Phys. , 86, 3471, (1987) ; J. Schnitker, K. Motakabbir, and P.
J. Rossky, Phys. Rev., 60, 456, (1988) .
[18] S. Huzinaga, Can. J. Chem., 73, 619, (1995) ; P. J. Hay and W. R. Wadt, J. Chem.
Phys., 82, 270, 284, 299, (1985) ; W. J. Stevens, H. Basch, and M. Krauss, J. Chem.
Phys., 81, 6026, (1984) ; W. J. Stevens, M. Krauss, H. Basch, and P. G. Jasien, Can.
J. Chem., 70, 612, (1991) など[19] K. Ohta, Y. Yoshioka, K. Morokuma and K. Kitaura, Chem. Phys. Lett. , 101,
12,(1983).
[20] S. Huzinaga and A. A. Cantu, J. Chem. Phys. , 55, 5543, (1971) ; Y. Sakai and S.
Huzinaga, J. Chem. Phys. , 76, 2537, (1982) ; S. Katsuki, Chem. Phys. Lett. , 127,
231, (1986) ; S. Katsuki and S. Huzinaga, Chem. Phys. Lett. , 152, 203, (1988) ; S.
Katsuki, Can. J. Chem. , 70, 285, (1992) ; S. Katsuki, J. Chem. Phys. , 98, 496,
(1993) ; S. Katsuki, Can. J. Chem. , 73, 696, (1995).
[21] N. Vaidehi, T. A. Wesolowski and A. Warshel, J. Chem. Phys. , 97, 4264, (1992) ; T.
A. Wesolowski and A. Warshel, J. Phys. Chem. , 97, 8050, (1993) ; T. A. Wesolowski
and A. Warshel, J. Phys. Chem. , 98, 5183, (1994).
[22] K. Honda and K. Kitaura, Chem. Phys. Lett. , 140, 53, (1987) ; K. Honda and K.
Kato, Chem. Phys. Lett. , 229, 65, (1994).
[23] J. N. Murrell, M. Randic and D. R. Williams, Proc. Roy. Soc. A, 284, 566, (1964) ;D. R. Williams, L. J. Schaad and J. N. Murrell, J. Chem. Phys. , 47, 4916, (1967).
[24] C. Muguruma, N. Koga K. Kitaura and K. Morokuma, Chem. Phys. Lett. , 224, 139,
(1994) ; C. Muguruma, N. Koga K. Kitaura and K. Morokuma, J. Chem. Phys. ,103, 9274, (1995).
[25] 例えば、W. L. Jorgensen, J. Am. Chem. Soc. , 101, 2011,(1979) ; W. L. Jorgensen,
J. Am. Chem. Soc. , 103, 335, (1981) ; W. L. Jorgensen, J. Chem. Phys. , 77, 4156,(1982) ; W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M.
L. Klein, J. Chem. Phys. , 79, 926, (1983) ; W. L. Jorgensen, J. D. Madura and
C. J. Swenson, J. Am. Chem. Soc. , 106, 6638, (1984) ; W. L. Jorgensen and C. J.
Swenson, J. Am. Chem. Soc. , 107, 569, (1985) ; W. L. Jorgensen and J. TiradoRives, J. Am. Chem. Soc. , 110, 1657, (1988) ; W. L. Jorgensen, J. M. Briggs and
M. L. Contreras, J. Phys. Chem. , 94, 1683, (1990).
43

第 6. まとめと展望

第四分科会・佐藤啓文[26] S. F. Boys and F. Bernardi, Mol. Phys., 19, 553, (1970)[27] J. H. van Lenthe, J. G. C. M. van Duijneveldt - van de Rijdt, F. B. van Duijneveldt,
Ab initio Methods in Quantum Chemistry II, 521 Ed. K. P. Lawley, John Wiley &Sons Ltd., 1987
[28] D. Rinaldi and J. L. Rivail, Theor. Chem. Acta, 32, 57, (1973) ; J. L. Rivail and D.
Rinaldi, Chem. Phys. , 18, 233, (1976).
[29] D. Rinaldi, M. F. Ruis-Lopez and J. L. Rivail, J. Chem. Phys. , 78, 834, (1983) ; D.
Rinaldi, B. J. Costa Cabral and J. L. Rivail, Chem. Phys. Lett. , 125, 495, (1986) ;D. Rinaldi, J. L. Rivail and N. Rguini, J. Comp. Chem. , 13, 675, (1992).
[30] J. H. McCreery, R. E. Christoﬀersen and G. G. Hall, J. Am. Chem. Soc. , 98, 7191,(1976).
[31] M. W. Wong, M. J. Frish and K. B. Wiberg, J. Am. Chem. Soc. , 113, 4776, (1991) ;M. W. Wong, K. B. Wiberg and M. J. Frish, J. Am. Chem. Soc. , 114, 523, (1992) ;M. W. Wong, K. B. Wiberg and M. J. Frish, J. Am. Chem. Soc. , 114, 1645, (1992); M. W. Wong, K. B. Wiberg and M. J. Frish, J. Chem. Phys. , 95, 8991, (1991).
[32] Y. Jeanvoine, F. Bohr and M. F. Ruiz-López, Can. J. Chem. , 73, 710, (1995). ; R.
J. Hall, M. M. Davidson, N. A. Burton and I. H. Hillier, J. Phys. Chem. , 99, 921,
(1995).
[33] S. C. Tucker and D. G. Truhlar, Chem. Phys. Lett. , 157, 164, (1989) ; C. J. Cramer
and D. G. Truhlar, J. Am. Chem. Soc. , 113, 8305, (1991) ; D. J. Giesen, J. W.
Storer, C. j. Cramer and D. G. Truhlar, J. Am. Chem. Soc. , 117, 1057, (1995).
[34] S. Miertuš, E. Scrocco and J. Tomasi, Chem. Phys. , 55, 117, (1981) ; S. Miertuš
and J. Tomasi, Chem. Phys. , 65, 239, (1982) ; G. Algona, C. Ghio, J. Igual and
J. Tomasi, J. Am. Chem. Soc. , 111, 3417, (1989) ; F. J. Olivares del Valle and J.
Tomasi, Chem. Phys. , 150, 139, (1991) ; M. A. Aguilar, F. J. Olivares del Valle andJ. Tomasi, J. Chem. Phys. , 98, 7375, (1993) ; M. A. Aguilar, R. Bianco, S. Miertus,
M. Persico and J. Tomasi, Chem. Phys. 174, 397, (1993) ; R. Cammi and J. Tomasi,
J. Chem. Phys. , 100, 7495, (1994) ; R. Cammi and J. Tomasi, J. Chem. Phys. ,101, 3888, (1994).
[35] H. L. Friedman, Mol. Phys. , 29, 1533, (1974).
[36] G. Karlström, J. Phys. Chem. , 92, 1315, (1988) ; G. Karlström, J. Phys. Chem. ,92, 1318, (1988) ; G. Karlström, J. Phys. Chem. , 93, 4952, (1989). M. Andersson
and G. Karlström, J. Phys. Chem. , 89, 4957, (1985)[37] H. Sato, and S. Kato, J. Mol. Struct. (THEOCHEM), 310, 67, (1994).
44

第 6. まとめと展望

第四分科会・佐藤啓文[38] A. Klamt and G. Schüürmann, J. Chem. Soc. , Perkin Trans. II, 799, (1993).
[39] T. N. Truong and E. V. Stefanovich, Chem. Phys. Lett. , 240, 253, (1995) ; E. V.
Stefanovich and T. N. Truong, J. Chem. Phys., 105, 2961, (1996).
[40] H. Hoshi, M. Sakurai, Y. Inoue and R. Chujo, J. Chem. Phys. , 87,1107,(1987).
[41] W. L. Jorgensen, Acc. Chem. Res. , 22, 184, (1989) ; B. J. Gertner, R. M. Whitnell,
K. R. Wilson and J. T. Hynes, J. Am. Chem. Soc. , 113, 74, (1991) ; S. E. Huston,
P. J. Rossky and D. A. Zichi, J. Am. Chem. Soc. , 111, 5680, (1989).
[42] S. Kato, and Y. Amatatsu, J. Chem. Phys., 92, 7241, (1990).
[43] A. Warshel, Computer Modeling of Chemical Reactions in Enzymes and Solutions,John Wiley & Sons, Inc., (1991).
[44] J. Gao and X. Xia, J. Am. Chem. Soc. , 115, 9667, (1993) ; P. Bala, B. Lesyng and
J. A. McCammon, Chem. Phys. , 180, 271, (1994) ; M. A. Thompson and G. K.
Schenter, J. Phys. Chem. , 99, 6374, (1995).
[45] S. Ten-no, F. Hirata and S. Kato, Chem. Phys. Lett. , 214, 391, (1993) ; S. Ten-no,
F. Hirata and S. Kato, J. Chem. Phys. , 100, 7443, (1994) ; M. Kawata, S. Ten-no, S.
Kato, and F. Hirata, Chem. Phys. Lett. , 240, 199, (1995) ; M. Kawata, S. Ten-no, S.
Kato, and F. Hirata, J. Am. Chem. Soc. , 117,1638, (1995). M. Kawata, S. Ten-no,
S. Kato, and F. Hirata, J. Phys. Chem., 100, 1111, (1996) ; M. Kawata, S. Ten-no,
S. Kato, and F. Hirata, Chem. Phys. , 203, 53, (1996).
[46] H. Sato, F. Hirata, and S. Kato, J. Chem. Phy., 105, 1546, (1996).
[47] T. Ishida, F. Hirata, H. Sato, and S. Kato, J. Phys. Chem., B, 102, 2045, (1998).

45

付 録A

汎関数微分とその分布関数理論への応用1２体分布関数 ρ(2) (r, r0 ) は２粒子間相互作用 v(r) に依存する汎関数であるのみならず、rと r0 の関数でもあるので、ρ(2) (r, r0 ) = ρ(2) [v; r, r0 ] と書かれる。この例からわかるように、
統計力学で扱われる多くの量はある関数の汎関数である。ここでは汎関数微分の基本的な事柄についてまとめる。

A.1

汎関数微分

Γ がある関数 u(x) の汎関数であるとし、それを Γ = Γ[u] と書くことにする。関数 u が
u + δu に変化した時、Γ が Γ + δΓ に変化したとすると、その変化分は δu の一次までとると次式の様になる。
δΓ = Γ[u + δu] − Γ[u]∫=dxA[u; x]δu(x)(A.1)

最後の式で汎関数 Γ[u] の関数 u(x) に関する汎関数微分、

A[u; x] ≡

δΓ
δu(x)

(A.2)

を定義した。２階の汎関数積分は次のようにして定義される。
∫

δA[u; x]δu(x0 )
δu(x0 )∫
δ2Γ
δu(x)δu(x0 )≡dx0
δu(x)δu(x0 )δA[u; x] =

dx0

(A.3)

より高階の汎関数微分も同様にして定義される。これらを利用すると、以下の汎関数テイラー展開を得ることができる。
∫

Γ[u] = Γ[u0 ] +
1
+
2!

∫∫

{

δΓ
dr
δu(r)
{

0

drdr

}

[u(r) − u0 (r)]
u=u0

δ2Γ
δu(r)δu(r0 )}[u(r) − u0 (r)][u(r0 ) − u0 (r0 )]
u=u0

+···

(A.4)
46

付 録 A. 汎関数微分とその分布関数理論への応用2

第四分科会・佐藤啓文

以下に便利な公式をまとめておく。

δ(cΓ)
δΓ
= c
δf (r)
δf (r)
δ(Γ1 Γ2 )
δΓ1
δΓ2
= Γ2
+ Γ1
δf (r)
δf (r)
δf (r)
( )
1
1 δΓ
δ
= − 2
δf (r) ΓΓ δf (r)δΓ1 (r)=δΓ2 (r0 )∫δΓ1 (r0 , r00 )=δΓ2 (r1 , r2 )
dr00

{

{

(A.5)
(A.6)
(A.7)

δ(r − r0 ) (if Γ1 = Γ2 )0
(if Γ1 is not a functional of Γ2 )(A.8)

δ(r0 − r1 )δ(r00 − r2 ) (if Γ1 = Γ2 )(A.9)
0
(if Γ1 is not a functionl of Γ2 )

δΓ1 (r) δΓ2 (r00 )= δ(r − r0 )
δΓ2 (r00 ) δΓ1 (r0 )(鎖則 (Chain Rule))(A.10)

ここで次節への準備となる具体例を挙げておこう。以下の汎関数を考える。
(∫

Γ[u] =∫=)Ndxk(x)u(x)···∫dx1 · · · dxN k(x1 ) · · · k(xN )u(x1 ) · · · u(xN )≡ (G[u])N(A.11)

この汎関数を一階汎関数微分すると、

δΓ
δΓ δG
=
δu(x)
δG δu(x)(∫= N k(x))N −1dxk(x)u(x)(A.12)となり、また同様にして２階微分はδ2Γ= N (N − 1)k(x)k(x0 )δu(x)δu(x0 )

(∫

)N −2

dxk(x)u(x)(A.13)となる。

A.2

分布関数理論への応用

以下の式で与えられる相互作用ポテンシャルを持つ系を考える。

VN (1, · · · , N ) =N∑

φ(i) +

N∑∑v(i, j)(A.14)i<ji=1

この式において、粒子の位置 qi を表すのに簡単のため i を用いている。v(i, j) は２粒子間相
互作用、φ(i) は外場を表す。この時系のグランドカノニカル分配関数は次式で与えられる。

Ξ =

∞
∑
1 ∫
N =0 N !

···

∫ ∏
N
i=1

z ∗ (i)N∏e(i, j)d1 · · · dN

i<j

47

(A.15)

付 録 A. 汎関数微分とその分布関数理論への応用3

第四分科会・佐藤啓文

ここで、z ∗ (i) = z exp[−βφ(i)](A.16)e(i, j) = exp[−βv(i, j)](A.17)である。ここで n− 体分布をもう一度書き下すと、

ρ

(n)

∫
∫ ∏
∞
N
N
∏
1 ∑
1
(1, · · · , n) =···z ∗ (i) e(i, j)d(n + 1) · · · dN
Ξ N =n (N − n)!
i=1
i<j(A.18)である。この式は式 (A.12) を参考にすると以下のように変形できる。

δΞ=δz ∗ (1)=∫∞∑N
N =1 N !

···

∫ ∏
N

z ∗ (i)

i=2

∞
1
1 ∑
z ∗ (1) N =1 (N − 1)!

N∏e(i, j)d(n + 1) · · · dN

i<j

∫

···

∫ ∏
N
i=1

z ∗ (i)N∏e(i, j)d2 · · · dN

(A.19)

i=1

これにより、z ∗ (1) δΞΞ δz ∗ (1)δ log Ξ=δ log z ∗ (1)ρ(1) (1) =(A.20)となる。一般的に次のように書けることを示すことができる。

ρ(n) (1, · · · , n) =1 ∗δnΞz (1) · · · z ∗ (n) ∗Ξδz (1) · · · δz ∗ (n)48(A.21)

付 録B
B.1

RISM 方程式の導出

局所密度場と相関関数

体積 V の中に N 個の分子からなる体系があるとする。この容器内の一点 r における α種の原子（相互作用点）の局所密度は Dirac のデルタ関数を使って次式で表現できる。

ν̂α (r) =N∑δ(r − rαi )

(B.1)

i=1

ρ̂α (r) の統計平均、
ρα (r) = hρ̂α (r)i(B.2)

は位置 r に α-種の原子を見いだす確立に比例し、１体（原子）分布関数と呼ばれる。もし体系が一様な流体であれば ρα (r) は α 原子の平均密度に一致する。これらの量を使って密度−密度相関関数を次式で定義する。

χαγ (r, r0 ) ≡ h(ρ̂α (r) − ρα (r))(ρ̂γ (r) − ργ (r0 ))i(B.3)

式 (B.3) に式 (B.1) 、(B.2) を代入し、和を i = j （分子内）と i 6= j （分子間）の項に分けると、式 (B.3) は次のように書き換えられる。
000χαγ (r, r0 ) = χ(0)αγ (r, r ) + ρα (r)ργ (r )hαγ (r, r )

(B.4)

0
0
ここで、χ(0)αγ (r, r ) および hαγ (r, r ) はそれぞれ分子内密度相関関数および分子間２体相関関数であり、次式で与えられる。
0
0
0
χ(0)
αγ (r, r ) = ρα (r)δ(r − r )δαγ + ρα (r)(1 − δαγ )sαγ (r, r )ρα (r)ργ (r0 )hαγ (r, r0 ) =

*
∑∑
i

+

δ(r − rαi )δ(r0 − rγj ) − ρα (r)ργ (r0 )

(B.5)
(B.6)

j6=i

ただし、sαγ (r, r0 ) は
0ρα (r)sαγ (r, r ) =*∑+δ(r − rαi )δ(r0 − rγi )

(B.7)

i

で定義される分子内相関関数であり、剛体分子の場合サイト α、γ 間の距離を Lαγ として

sαγ (r, r0 ) =1
δ(|r − r0 | − Lαγ )
2
4πLαγ

(B.8)

の表現を得る。

49

付 録 B. RISM 方程式の導出

第四分科会・佐藤啓文

RISM(SSOZ) 方程式

B.2

体系に場 φα がかかっている時、そのポテンシャルエネルギーを次式で定義する。

Φext =∑∫

drρα (r)φα (r)

(B.9)

α

この場に対する ρα (r) の応答は容易に求めることができて、次式を得る。

δρα (r)= χαγ (r, r0 )
δ[−βφγ (r0 )](B.10)

式 (B.10) の逆関係により直接相関関数 c(r, r0 ) を次式で定義する。

δ[−βφα (r)]= χαγ (r, r0 )−1δργ (r0 )1=
δαγ δ(r − r0 ) − c̄αγ (r − r0 )
ρα (r)(B.11)

0
同様に、分子内の直接相関関数 c(0)αγ (r, r ) を次式で定義する。
0 −1
χ(0)
=
αγ (r, r )10δαγ δ(r − r0 ) − c(0)αγ (r − r )ρα (r)(B.12)

式 (B.11) 、(B.12) より、

δ[−βφα (r)]0 −10= χ(0)− {c̄αγ (r − r0 ) − c(0)αγ (r, r )αγ (r, r )}0δργ (r )(B.13)

式 (B.11) および式 (B.13) を汎関数微分の鎖則（式 (A.10) を参照）∑∫ηδρα (r) δ[−βφη (r00 )] 00
dr = δαγ δ(r − r0 )
δ[−βφη (r00 )] δργ (r0 )(B.14)に代入し、式 (B.4) 、(B.5) 、(B.6) を考慮すると φα (r)（一様な流体）の場合およびに関する次の積分方程式を得る。

ρhρ = ω ∗ c ∗ ω + ω ∗ c ∗ ρhρ(B.15)

ただし、∗ は行列の乗算と畳み込み積分を表し、ρ は平均密度の行列である。h、c、ω は原
子対を添え字とする行列であり、それらの要素はそれぞれ式 (B.6) および以下の式によって与えられる。

cαγ (r) ≡ c̄αγ (r) − c(0)αγ (r)(B.16)

ωαγ (r) ≡(B.17)χ(0)αγ (r)

ただし、r は α − γ サイト間の距離を表す。式 (B.15) を RISM 方程式、あるいは SSOZ 方
程式と呼ぶ。

50