測定法 Q & A 電気化学そこを知りたい，議論したい
量子化学計算による酸化還元電位の算出法
Redox Potential by First-principles Calculation
測定法小委員会

Modified on December 13, 2003
Ｑ 30. 「量子化学計算による酸化還元電位の算出法特に，（簡単の為に，以後 z = +1, n = 1 とします。）の酸化還元
種々の量子化学計算手法，電気化学系に応用した場合の問題 電位を左側のセルの水素電極を基準にとって，開回路電圧点等について，初心者にとっても分かり易いように説明して (OCV) すなわち平衡電位として計算します。Cu| Pt | H2ほしい。」(P = P ◦ = 105 Pa ) | H+ (aH+ = 1) || Ox+ / Red | Pt’ |Cu’ のセルでの開回路電圧 E は，これまでこの Q&A シリー
ズで詳細に検討されているように [2, 3], 半電池反応 Eq.(1)およびA30:酸化還元電位を理論的に求めることは，酸化還元電1* H2(2)
H+ + e− )位が何により決定されているのかという基本的な問題を明2らかにする点だけにとどまらず，例えば二次電池の起電力を理論的に予測するという実用的な観点からも興味深いよう の半電池反応に対する平衡関係からに思います。
酸化還元電位の算出法，その近似の程度により幾つかの方法があるかと思います。狭義な意味での量子化学 (QC) 計算（すなわち，LCAO(Linear Combination of Atomic Orbital)基底系でのクラスター計算）ではなく，平面波基底系を用いた周期格子系に対するバンド計算 (BC: Band Calculation)も含めた第一原理計算（第一原理計算とは量子力学の基礎方程式から出発すると言う意味で，ab initio 計算とも言います。狭い意味で使う場合は，密度汎関数理論に基づく電子状態計算のことです。) でいかに計算するかについて最近の研究も含めて解説します。

CuEeφPt

ここで，φi は i でのポテンシャル，µ◦i は i の標準状態での化学ポテンシャル，kB は Boltzmann 定数，T は温度，ai はi の活量です。化学ポテンシャルはモル数単位ではなく粒子数単位で示されています [4]。平衡電位なので溶液内の IRドロップはゼロで，もし液間電位が無視できるとするなら，
φH+ = φOx+ であり，また µ◦Pt/e = µ◦Pt0 /e ですのでEPt e-H+1/2 H 2ReOx+(3)1kB T− µ◦H2 −ln(PH2 /P ◦ )(4)22= µ◦Ox+ + kB T ln aOx+ + eφOx+ + µ◦Pt0 /e−µ◦Re − kB T ln aRe(5)

eφPt0Cu’V= φCu0 − φCu = φPt0 − φPt= µ◦H+ + kB T ln aH+ + eφH+ + µ◦Pt/ee-=1 ◦1(µ + − µ◦Re − µ◦H+ + µ◦H2 )e Ox
2kB TaOx+ (PH2 /P ◦ )1/2+lneaRe aH+(6)標準電極電位 E ◦ は，すべての活量が１ mol dm−3 で水素圧が 105 Pa である時の E として定義されるのでPt’E◦ =H2 gas

11 ◦(µ + − µ◦Re − µ◦H+ + µ◦H2 )e Ox2(7)となります。化学ポテンシャルの定義から, 標準状態での化学ポテンシャルは，µ¶∂G◦µi =(8)∂Ni ◦Figure 1: Schematic diagram of the cell to measure thestandard electrode potential of the redox reaction.

まず，計算方法を述べる前にその定義を復習しましょう。 標準状態で i を無限小増減させた時（すなわち，１粒子増減
このシリーズの Q7 の標準電極電位測定法 [1] にあるように させた時）の自由エネルギー変化です。水素ガスの標準状態
Fig.1 に示したような電極で右側のセルでの酸化還元反応
での，µ◦H2 は，剛体回転子-調和振動子モデルで統計力学で求められているのでそれを用い [5]，水素分子の電子エネル
Oxz + ne− *(1) ギーは経験的なパラメータのない第一原理計算から求めま) Rez−n1す。電子エネルギーの基準の取り方は任意であるので，例え 原理計算でこれらの方法により溶媒和自由エネルギーを計算ば，電子を持たない真空中にあるプロトンをエネルギーゼ した例はないようです。Li+ イオンの場合 Lyubartsev らはロの基準にとります。それぞれの化学種の標準状態における [22]，第一原理分子動力学計算の結果からリチウムと水の間
化学ポテンシャル µ◦Ox+ , µ◦Re , µ◦H+ を理論的に求めることは， の有効ポテンシャルを求め，それを用いて 1000 水分子のよ第一原理計算に基づいて系の自由エネルギーを高精度に求 り大きな系で熱力学的な積分をすることで溶媒和自由エネルめることであり，後述するように困難な問題です。扱う系に ギーを求めました。溶媒和自由エネルギーは，-419 kJ/mol
よって近似のレベルも異なってくるので以下は個々の場合を となり，これまで報告された理論計算に比較して，より実験考えます。
結果-480 kJ/mol[8] に近い値となりました。
また，溶媒和のエントロピー項（プロトンで+30kJ /酸化還元電位がリチウムイオン電池（酸化物とグラファイト）のような固相反応で決定される場合は，リチウムイ mol ) を考慮せずにエンタルピー (プロトンで-1094 kJ / mol)
オン濃度 (サイトの占有率）θ での全エネルギー U (θ)（即ち だけを議論するときは，Fig.2 に示したよう各状態間の内部内部エネルギー。固体の場合は内部エネルギーはエンタル エネルギー差で議論できるかもしれません。この場合，(7)ピーにほぼ等しい。）を θ について何点か計算し, その微分 式での化学ポテンシャルは系のエンタルピーをイオンの数で値 dU/dθ より化学ポテンシャルを求めます。97 年に報告さ 割ったものにほぼ等しいと近似するということです。即ち，れた Aydinol らの計算では，Li 基準の 0 ≥ θ ≥ 1 における 水素分子・水素原子の全エネルギー，水素イオン (又は酸化平均の電位を U (LiMO2 ) - U (MO2 ) - U (Li) で計算しまし 体・還元体) ＋電荷補償用の background 負電荷系の全エネた [6]。(ここで, M はコバルト等の金属を意味する。）この ルギー（計算には周期境界条件を用いるので，基本セル内で際，簡単のために U は絶対零度の値を用い，エントロピー 電荷が中性でないとよほど大きなセルを用いない限り全エ項は熱エネルギー程度（30 meV) なので無視したようです ネルギーが収束しなくなります。従って，イオンの電荷に等が，実験で求めた開回路電圧 (OCV) と一致しました。また， しい逆符号の均一な電荷を基本セル内に人為的に導入してＭの種類を変えて OCV を理論的に予測しそれを実験的にも 計算を行ってます。他の方法では，同数のカウンターイオン確認しました [7]。このときの計算法はバンド計算 (BC) タ を系に導入する方法もあります。溶質が電荷を持たない場合イプの計算です。BC タイプの計算法についての詳細は後ほ は必要ありません。）, 液体の水の全エネルギー，水素イオン（又は酸化体・還元体）＋電荷補償用の background 負電ど述べます。
◦
難しいのは液相での酸化還元反応です。 µi は活量が 1 の 荷＋水系の全エネルギーが計算できたとします。（これは推時（活量係数が 1 であれば濃度が 1 mol dm−3 ）の液相での 測ですが，上で紹介した計算において系のエネルギーを議論化学ポテンシャルを求める必要があります。この場合，溶媒 している報告例はありません。系の水分子数が少ないので溶相が液体であるには温度は融点以上でなければならず，絶対 媒和エネルギー（エンタルピー）は収束してないかもしれま
零度の ab-initio 計算では不十分です。また，真空中でイオ せん。）水素の場合と今着目している酸化体/還元体の場合ンの周りに溶媒分子を配置するようなクラスター量子化学 でそれぞれこれらの系の全エネルギーを計算しその差をと計算が数多く報告されてますが，溶媒分子の数を増加して れば，水素基準の酸化還元電位を求めることが可能なはずいった時，エネルギーが収束しないので，溶媒和エネルギー です。
現実に，酸化還元電位を求めるには，より粗い近似，すなを量子化学計算で理論的に求めることは困難であることがわち溶媒の効果を全く無視する(A) か溶媒を連続体で近似
以前から知られてました [8]。
するか（Ｂ），量子系と古典系を組み合わせる (C) 方法がと溶液中での溶媒和の理論計算では，実験結果に fit させたられます。
経験的なパラメータを用いた分子動力学計算により行われ(A) では真空中に酸化体又は還元体を置きi)HOMOてきており，fitting された環境以外へのパラメータの移植
LUMO レベルの差が化学ポテンシャルの差になると近似し可能性（transferability）が問題となり計算の信頼性に乏し酸化還元電位を見積もる。
ハートリーフォック法
いことが度々指摘されてきました。現在の最先端の第一原（HF 近似で，電子間の相関を考慮して物性値を計算をする
理計算では，問題となっていた原子・分子間力を経験的なときには さらに Møller-Plesset の摂動計算等が必要です。
パラメータのない第一原理計算から求め，系の統計的な熱密度汎関数理論では，電子間の相互作用に相関の項も含め平均量を分子動力学計算から導くことが可能になりました。て計算しているので摂動計算は必要ありません。）
では，その詳しい手法は，Marx と Hutter の 150 ページにおよぶ計算された軌道の１電子エネルギーはその軌道から電子をreview に詳しく [9]，また，これまで計算された結果は Tse取り去るのに必要なエネルギー（Koopmans の定理）に対の review[10] にもありますのでそれらをご参照下さい。報告
応するので HOMO-LUMO エネルギーレベルの計算から直された溶媒和の第一原理分子動力学計算には以下のものがあ接求めることができます。
ります。水 [11, 14]，溶媒和電子 [27]，プロトン・水酸化物ii) 密度汎関数法 (DFT) では，１電子の軌道エネルギー++
イオン [12, 13, 26]，水の解離定数 [21]，Li [22], Na [20],は系から無限小の電荷を取り去った時の全エネルギーの変++−2+2+
K [19], NH4 [18], Br [24], Be [16], Ca [23], グルコー化に相当します。(Janak の定理）DFT では，むしろ，系のス [17] の溶媒和，強酸 HCl[15], HF[25] の溶媒和などです。
全エネルギーが精度よく計算できるので酸化体の全エネル
それらの計算では，32-120 分子を周期境界条件をもつ基本ギーと還元体の全エネルギーの差を標準化学ポテンシャルのセル内に置き (水分子 56 個に対しイオン 1 分子でほぼ 1 mol差に近似して酸化還元電位を求めます。
−3dm になります。），単位時間ステップあたり 0.15fs で 10psi),ii) の両方の手法を用いて宇恵氏らは最近注目されていまでの計算を行ってます。イオン周りの水（酸素）の動径分る常温溶融塩の酸化電位を量子化学計算を用いて見積もり，布，配位数，配位角，溶媒和ダイナミックス，振動スペクト実験結果とのよい相関を示しました。[28]ルが実験結果とよく一致しました。ただし，溶媒和自由エネ(A) の計算法では，標準電極電位の絶対値は求められない
ルギーに関してはその計算が簡単ではありません。拡張アン
ことを，サンブル法 [22] や，始状態から終状態までの間のいくつかの中間状態における平均力ポテンシャルを計算してそれを積分する必要があります [10]。残念ながら，現在のところ第一

1Cl2 + e− *) Cl− ,22E ◦ = +1.3595V(9)その他，(i) 電極からのトンネリングによる電子移動，(ii)電気２重層効果等の電極表面での効果も考慮しなくてはならないでしょう。その意味で，電気化学を理論的に取り扱うには空間スケールでは電子レベルから連続体モデルまで，時間レベルでは fs から秒・分まで関わるマルチスケールでの理論解析が必要です。酸化還元電位が理論から正確に求めることが不可能ではなくなってきた現在，電気化学は Nernst以来のパラダイムシフトを迎えているのかも知れません。
計算ソフトウェア
密度汎関数理論 (DFT) に基づいた計算が主に行われています。Schrödinger 方程式に相当する Kohn-Sham 方程式を３次元で解くわけですが，実空間で解くかかフーリエ変換された逆空間で解くかで QC 系と BC 系に大きく分かれます。
自分で一から作成できれば良いのですが，数値計算の手法・プログラミングには過去の膨大な蓄積があるので以下に示すようなソフトウェアーを入手して使用するので望ましいと思います。可能な限り，ソースコード（基本的にはFortran90）が入手できるようなソウトウェアが良いように思います。
プログラムや数値計算の手法を完全に理解する必要はないかもしれませんが，第一原理計算でも種々の近似がなされておりその計算の原理およびその限界を十分理解した上でないと大きな間違いをする危険性もあります。
（たとえば，HF 法では一酸化炭素の双極子の向きは逆になります。もちろん，基底系の精度を上げ DFT を用いれば実験値に近づきます。）例えば，F. Jensen ”Introduction to ComputationalChemistry” (Wiley, 2000) では，（非理論家で）新たに計算を始めようとする者に対して良い解説を与えていますので，一読をお薦めします。
第一原理計算の精度は，1-10 meV (0.1 - 1 kJ mol−1 ) 程度だと思います。

Figure 2: Schematic diagram of the periodic unit cell tocalculate the solvation energy of ion.の例を挙げて示しましょう。計算は，ガウシャン 98 を用いB3LYP/6-311+G**レベルで行いました。
（記号の意味は，後述します。）全電子数が偶数の時，スピン多重度 2S+1 = 3 に，量子化学計算 (QC) 系：（LCAO 基底, クラスター計算）
奇数の時は 2S + 1 = 2 にしました。(1/2)H2 , H+ ，(1/2)Cl2 ,−１）ガウシャン(Gaussianinc.)[32]：最近 03 バージョン
Cl の原子あたりの真空中での内部エネルギーはそれぞれになり周期系の計算等かなり拡張されたようです。どのよ
eV 単位で，-15.912, 0.0, -12522.793, -12525.501 となりました。Eq.(7) より，E ◦ は，-13.204 V となり，実験結果とは うな計算をさせるのかを指定し（１行），電荷，多重度（１全く値が異なります。計算では，水素原子のイオン化エネ 行），原子座標（原子数行）を指定するだけです。初心者でルギー 13.664(実験値 13.598 eV)，塩素原子の電子親和力は も入力をグラフィックインターフェースを用いればすぐに計
3.723 eV（実験値 3.617 eV) とよく一致しているので，問 算可能で使いやすいように工夫されてますが，実際，中でな題は溶媒和のエネルギーにあります。実験で求めたプロト にを計算しているかは，一般のユーザーは闇のなかにある
ン及び塩化物イオンの溶媒和自由エネルギーは, -10.944 eV, と言って良いかもしれません。近年，ハイエンドパソコン用-3.596 eV [8] であるので，その効果を入れると標準電極電位 の CPU とワークステーション・スパコンの CPU の差はな
は+1.696 V となり実験値+1.3595 V に近い値となります。 くなりつつあります。その意味でガウシャン等の第一原理計以上のことから，第一原理計算において溶媒和のエネルギー 算もお手元の PC 上で十分計算可能です。対象とする物質にを見積もることが酸化還元電位の絶対値を見積もる上では もよりますが，筆者の研究室では B3LYP タイプ密度汎関数法を用い，基底関数はなるべく精度の高いもの 6-311+G**必要不可欠であることが明らかとなりました。
B) では，分子を誘電率をもつ連続媒体に挿入し量子化学 を用いています。
分子系に密度汎関数理論が適用されるようになったのは，計算を行うもので，反応場 (RF:reaction field) 自己無撞着汎関数に，局所密度近似（
LDA: Local Density Approxima反応場 (SCRF:self-consistent reaction field) モデルがありtion)に加えて一般化された勾配近似(GGA: Generalizedます。ガウシャンでもこの種の計算が可能です。[29]Gradient Approximation) が導入されて計算の精度が良くC) では，例えばイオンと第一水和圏にある分子のみを量なってから以降のことです。電子の交換・相関相互作用の汎子化学系で解き，その外側にある分子を古典的な分子動力関数の一つである B3LYP は，Becke’s 3 parameters Hybrid
学（又はモンテカルロ）計算で行う QM(Quantum mechanFunctional using Lee, Yang, and Parr correlation Functuics)/MM(Molecular Mechanics) 計算があります。量子系とional という意味で，Hartree-Fock(HF) の計算に含まれてい古典系をどのように結合させるかという困難な問題もありまる交換項と DFT での交換・相関項を混合して汎関数をつくっすが，実用的には今後この方法が採用されていくでしょう。
XCたもので，そのエネルギー EB3LYP
はイオンの水和自由エネルギー [30]，酵素の酸化還元電位 [31]XC
X
X
X
を QM/MM 法で求めた報告が最近なされています。
EB3LYP = c0 EHF + (1 − c0 )ESlater+ c1 ∆EBecke

3

X
ESlater=X∆EBecke=xσ=CC+c2 ELYP+ (1 − c2 )EVWNµ ¶1/3 Z
2 3−drρ4/3[33]3 4π4/3XZ
ρσ x2σ−βdr
1 + 6βxσ sinh−1 xσσ

ρ−4/3|∇ρσ |,σβ = 0.0042au

謝辞
この文章を書くにあたり御議論頂きました加納健司先生(京大院農)，前田耕治先生 (京都工繊大院物質)，垣内 隆先
生 (京大院工) に感謝致します。
キーワード酸化還元電位，学， 量子化学計算[34]

第一原理計算，分子動力h 文責 i 山本雅博 （京都大学大学院工学研究科）です。ここで，X は交換，C は相関を意味します。3 つのパラメータは，Becke[35] によって c0 = 0.20, c1 = 0.72, c2 = 0.81と与えられています。ρ は電子密度，σ は電子スピンの上向き，CC
文献の汎関数 E C [ρ(r), ∇ρ(r)], EVWN
下向きを示します。ELYPについては紙面の関係で省略しますが，興味のある方は参考
論文 [36, 37] と [38, 39] をそれぞれ参照下さい。
基底関数 6-311+G**は，次の意味を持ちます [40]。内 References
殻軌道に６個のガウス型関数，価電子には 3/1/1 の３重の[1] 北村房男， Electrochemistry, 70, 295 (2002).
大きさの異なる分割価電子軌道，**は d 分極関数とｐ分極[2] 加納健司，大堺利行，桑畑 進，Electrochemistry, 70関数を加え，さらに＋は大きく広がったｓまたは p 軌道を, 358 (2002).
加える事を意味します。計算で得られた結果を，基底系に分[3] 垣内 隆， Electrochemistry, 70 , 821 (2002).
解して議論するような特殊なことをするのでなければ，基底系が適用できる範囲（6-311+G**の場合水素から Kr まで）， [4] 通常はモル数で定義しますが，分子シミュレーションの場合，計算時間にほぼ３乗で比例する全クラスターでの基底関数
何個単位の分子数で議論することが多いので，粒子数当たりでの数（6-311+G**の場合炭素は 22 個の基底関数），基底関定義しました。アボガドロ数を NA とすれば，R = kB T, F =NA e となります。
数の種類に依存する計算の精度（実験値からの平均絶対偏差 (MAD: Mean Absolute Deviation from experiment) が[5] 例えば，D. A. McQuarrie, Statistical Mechanics (Univ. Sci.
Foresman とÆ. Frish の本に MAD 値として与えられていまBook) 2000 p.108
Eq.(6-56).
す。B3LYP/6-311+G**では 3.1 kcal / mol です。HF/STO[6] M. K. Aydinol, A. F. Kohan, G. Ceder, K. Cho, and J.
3G では，93.3 kcal / mol ですから 30 倍の精度差がありま
Joannopoulos, Phys. Rev.B 56, 1354 (1997).
す。）を基準にして，計算コストと計算精度の trade-off でDFT の関数系/基底系を選択されるのが良いかと思います [7] G. Ceder, Y. -M. Chiang, D. R. Sadoway, M. K. Aydinol,
Y. -I. Jang, and B. Huang, Nature, 392, 694(1998).
[29]。
[8] Y. Marcus, Ion Solvation, Wiley 1985.
２）その他 GAMESS (General Atomic and MolecularElectronic Structure System, M. S. Gordon research group, [9] D. Marx and J. Hutter, In Modern Methods and Algorithms
Iowa State Univ., Academic では free，ソースコード有)[41]of quantum Chemistry, Ed. J. Grotendorst, 1, pp.301-449(John von neumann Institute for Computing, Jülich), NIC
が QC 系のソフトウェアではあります。
Series vol. 1 2000.
[10] J. S. Tse, Annu. Rev. Phys, Chem., 53 249(2002).
バンド計算 (BC) 系:(平面波基底，周期境界セル）
バンド計算 (BC) 系ソフトウェアは，従来研究者の間で [11] K. Laasonen, M. Sprik, M. Parrinello, and R. Car, J.
opensource として出回っていたものをユーザーフレン
Chem. Phys. 99, 9080(1993).
ドリーな形に書き換えられたヨーロッパ系ソフトが使いや[12] M. Tuckerman, K. Laasonen, M. Sprik, and M. Parrinello,
すいようです。ウルトラソフト擬ポテンシャル法としては
J. Chem. Phys. 103, 150(1995).
VASP(Vienna Ab-initio Simulation Package) [42], 内殻レベルまで考慮した全電子計算である FLAPW(Full-potnetial [13] M. Tuckerman, K. Laasonen, M. Sprik, and M. Parrinello,
J. Phys. Chem. 99, 5749(1995).
Linearized Augumented Plane Wave) 法の計算コードではWIEN2k[43] があります。どのような計算を行っているか [14] M. Sprik, J. Hutter, and M. Parrinello, J. Chem. Phys.
105, 1142(1996).
を知るには，固体物理の知識が少々必要ですが，それぞれのマニュアルの参考文献から遡っていけば解ると思います。 [15] K. E. Laasonen and M. L. Klein, J. Phys. Chem. A 101,98(1997).
これらの計算を行っている研究者が近くにおればその方に相談頂くのが一番手っ取り早いように思います。最近では， [16] D. Marx, M. Sprik, and M. Parrinello, Chem. Phys. Lett.
自力で BC 系の計算を行い実験結果と詳細な比較をされて273, 360(1997).
いる電気化学系の研究室もあると伺ってます。[44] ab-initio[17] C. Molteni and M. Parrinello, J. Am. Chem. Soc. 120,
擬ポテンシャル法及びウルトラソフト擬ポテンシャル法の2168(1998).
解説は http://fm.ehcc.kyoto-u.ac.jp/japanese.html の selfstudy Notes にもあります。平面波基底系であるので，計算 [18] F. Bruge, M. Bernasconi, and M. Parrinello, J. Am. Chem.
Soc. 121, 10883(1999).
は LCAO 系基底と異なり平面波の数を増やしていけば単純に収束していきます。その他，計算に用いる逆格子点の数， [19] L. M. Ramaniah, M. Bernasconi, and M. Parrinello, J.
Chem. Phys. 111, 1587(1999).
FFT 格子点の数，計算するバンドの数，用いる密度汎関数法の種類 (LDA,GGA)，電子ポテンシャルの性質に精度は依 [20] J. A. White, E. Schwegler, G. Galli, and F. Gygi, J. Chem.
Phys. 113, 4668(2000).
存します。

4

[21] P. L. Geissler, C. Dellago, D. Chandler, J. Hutter, and M.
Parrinello, Science 291, 2121(2001).
[22] A. P. Lyubartsev, K. Laasonen, and A. Laaksonen, J.
Chem. Phys. 114, 3120(2001).
[23] I. Bako, J. Hutter, and G. Palinkas, J. Chem. Phys. 117,9838(2002).
[24] S. Raugei and M. L. Klein, J. Chem. Phys. 116, 196(2002).
[25] A. J. Sillanpaa, C. Simon, M. L. Klein, and K. Laasonen,
J. Phys. Chem. B 106, 11315(2002).
[26] M. E. Tuckerman, D. Marx, and M. Parrinello, Nature 417,925(2002).
[27] M. Boero, M. Parrinello, K. Terakura, T. Ikeshoji, and C.
C. Liew, Phys. Rev. Lett. 90, 226403 (2003).
[28] M. Ue, A. Murakami, and S. Nakamura, J. Electrochem.
Soc., 149 A1572(2002).
[29] J. B. Foresman and Æ. Frish, in Exploring Chemistrywith Electronic Structure Methods (2nd Ed.) (Gaussain Inc.
1996) Chap. 10.
[30] T. Shoeib, G. D. Ruggiero, K. W. M. Siu, A. C. Hopkinson,and I. H. Williams, J. Chem. Phys. 117, 2762(2002).
[31] M. H. M. Olsson, G. Y. Hong, and A. Warshel, J. Am.
Chem. Soc. 125, 5025(2003).
[32] G03:http://www.gaussian.com/ ,G98:http://www.lct.jussieu.fr/manuals/Programmes/Gaussian98/index.htm[33] J. C. Slater, Quantum Theory of Molecular and Solids.
Vol. 4: The Self-Consistent Field for Molecular and Solids(McGraw-Hill, New York, 1974).
[34] A. D. Becke, Phys. Rev. A, 38 3098 (1988).
[35] A. D. Becke, J. Chem. Phys., 98 5648 (1988).
[36] C. T. Lee, W. T. Yang, and R. G. Parr, Phys. Rev. B, 37785 (1988).
[37] B. Miehlich, A. Savin, H. Stoll, and H. Preuss, Chem. Phys.
Lett., 157 200 (1989).
[38] S. H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. 58
1200 (1980).
[39] R. M. Dreizler and E. K. U. Gross, ”Density FunctionalTheory” (Springer, Berlin 1990) Chap. 7.
[40] T. Clark, A Handbook of Computational Chemistry (Wiley,New York 1985) Chap. 5.
[41] http://www.msg.ameslab.gov/GAMESS/GAMESS.html[42] http://cms.mpi.univie.ac.at/vasp/[43] http://www.wien2k.at/index.html[44] 中山，内本，私信.

5