http://www.chem.ccu.edu.tw/~groupweb/
由胡維平教授所開的各門大學部及研究所課程可參考下列分頁:
為電子的動能,電子與原子核之間的位能,電子間的排斥位能,原子核的
動能,以及原子核之間的排斥位能。在 Born-Oppenheimer 假設下我們可以忽略倒數第
二項,而最後一項在固定原子核位置下是常數,通常以 VNN 來代表。因此,只包含電
子座標的薛丁格方程式可
1
量子化學計算
2009/10/14
背景
西元1929 年 P. A. M. Dirac 曾表示:『整個化學領域背後的物理定律之數學理論已
經全然知曉,困難的地方在於實際應用這些理論時所推衍出的方程式過於複雜而難以求
解。』的確,量子化學理論在1920 年代末期已然成型,但其在化學研究上的廣泛應用
卻要等到大約50 年後,這是因為量化理論中所需求解的薛丁格方程式太過複雜,除了
最簡單的分子外,在有限的計算資源下,對一般化學系統應用上需要做太大的簡化,所
以通常最多也只能提供定性的預測,而且在使用上需要非常專業的量化訓練及計算能
力,因此量化計算在80 年代以前只是非常少數理論物理化學家的研究工具。
80 年代開始,情況有了明顯的改變。由於超大型積體電路的發展,以往大型笨拙
而難以操作的大型電腦逐漸被當時所謂的迷你電腦及高速工作站 (Vax, Sun, Digital) 取
代,以往冷冰冰深奧的電腦作業系統也逐漸轉成功能強大但較為友善的 Unix 系統並且
開始提供圖形化的人機介面。 在理論方法上,許多準確度高且有效率的近似方法陸續
被開發出來,在計算軟體上,由John Pople 在70 年代所開發出的 Gaussian 程式已經
包含非常豐富的功能,並被廣泛應用在化學研究上,其他如 AMPAC, MOPAC 等程式
包含各種非常有用的semi-empirical 方法也廣受化學家的歡迎。
到了90 年代末期,個人電腦的速度大幅的改進 (Pentium II,III),性能上與傳統的
工作站的差距已經非常有限,在作業系統上功能完整的 Window 2000 及 Linux 使得個
人電腦開始成為量化計算的利器。同時,許多功能強大容易使用的量化計算的軟體也陸
續支援 Window 及 Linux,如 Gaussian, Spartan, Q-Chem, Hyperchem 等。
今天,以台幣約四-五萬元就可以買到一套功能極為強大的個人電腦工作站,包含
最新四核心技術的處理器,8GB 高速記憶體,以及 600 GB 以上的硬碟空間,配合
Window XP (or Vista) 或基本上是免費的 64 位元Linux 作業系統,其運算效能遠超過
十年前需要上百萬美金的超級電腦。因此,由於近年來高速運算資源的普及以及有效率
的理論方法的發展,計算化學已經逐漸成為各種領域化學研究當中不可或缺的一部份,
相信在短期內也將成為是當代主修化學相關領域的學生的基本訓練之一。
2
量化計算理論
對任意分子系統而言其總能量 operator (Hamiltonian) 可寫成
Σ Σ Σ Σ Σ
> >
= − ∇ −+ −∇ +
A B A B
A B
A
A
i A i A i j i j A
A
i
i
e R
e Z Z
r m
e
r
e Z
m
H
0 ,
2
2
2
0 ,
2
0 , ,
2
2
2
2 4π
1
2 4π 4π
ˆ
ε ε ε
h h (1)
等號右邊分別為電子的動能,電子與原子核之間的位能,電子間的排斥位能,原子核的
動能,以及原子核之間的排斥位能。在 Born-Oppenheimer 假設下我們可以忽略倒數第
二項,而最後一項在固定原子核位置下是常數,通常以 VNN 來代表。因此,只包含電
子座標的薛丁格方程式可寫成
total ele NN
ele ele ele ele
E E V
H ψE ψ= +
ˆ=
(2)
其中 ele Hˆ為 (1) 式中的前三項,Etotal 又稱為 Born-Oppenheimer energy。Eq. (2) 對一
般化學分子而言還是太複雜,無法求得解析解,因此我們必須對薛丁格方程式做一些簡
化或提出一些額外的假設來求得 Eq. (2) 的近似解。
(1) Hartree-Fock 理論
通常近代量化計算的出發點是 Hartree-Fock 理論,也就是用一個行列式 (slater
determinant) 來代表波函數。對於一個含有 2n 個電子的 closed-shell 系統,slater
determinant 可寫成
(2 ) (2 ) (2 ) (2 ) (2 ) (2 )
(2) (2) (2) (2) (2) (2)
(1) (1) (1) (1) (1) (1)
(2 )!
1
1 1 2 2
1 1 2 2
1 1 2 2
n n n n n n
n
n n
n n
n n
slater
φ α φ β φ α φ β φ α φ β
φ α φ β φ α φ β φ α φ β
φ α φ β φ α φ β φ α φ β
LL
M
M
LLL
LLL
Ψ = (3)
其中 φ1…φn 為填滿之正交分子軌域 (orthogonal molecular orbitals)。為什麼要這麼寫
呢?首先,我們假設在多電子分子中電子還是可以看成在分子軌域中獨立運動,而任一
個電子其所受的有效位能為分子內所有的原子核以及其他電子的平均庫倫作用力;因
3
此,近似波函數本來可寫成所有自旋軌域的乘積,但電子是無法區分 (indistinguishable)
的基本粒子,而且苞利定律要求電子的波函數在經過任二個粒子座標對調後必須變號
(antisymmetric),因此一種最簡單且具有一般性的作法就是將波函數寫成一個如Eq.(3)
之slater determinant,基本上就是將波函數寫成各種自旋軌域乘積的 antisymmetric 線
性組合,2n 個電子在不同的軌域中可有任意的排列,因此共有 (2n)! 項。由近似波函
數 (3) 所得的能量期望值為:
∫ Σ Σ Σ( )
= = =
= Ψ Ψ = + −n
i
n
j
ij ij
n
i
Eslater slater Hele slater d Hii J K
1 1 1
* ˆτ 2 core 2 (4)
1
0 1,
2
1
2
core (1)
2 4π(1) * φ τ
ε
φ d
r
Z
m
H i
A A
A
e
i ii ⎟⎟⎠⎞⎜⎜⎝⎛= ∫ −∇ −Σ h (5)
1 2
0 12
2
(1) (2)
4π(1) * (2) * φ φ τ τ
ε
φ φ d d
r
J e j i j i ij ⎟⎟⎠⎞⎜⎜⎝⎛= ∫ (6)
1 2
0 12
2
(1) (2)
4π(1) * (2) * φ φ τ τ
ε
φ φ d d
r
K e i j j i ij ⎟⎟⎠⎞⎜⎜⎝⎛= ∫ (7)
其中 Hcore 是代表個別電子的動能以及它與所有原子核間的作用力,如果電子間沒有
作用力,則 (4) 式中的Hcore 項就是電子的總能量, r1,A 代表電子與原子核A 間的距
離。式 (4), (6) 中的 J 代表在二個分子軌域間的傳統電子-電子庫倫作用力,也稱為
庫倫積分 (Coulomb integrals) 或庫倫能量。式 (4), (7) 中的 K 雖然類似 J 但並沒有古
典力學的相對值,它的來源是苞利定律中要求波函數必須為antisymmetric,K 也稱為交
換積分 (Exchange integrals) 或交換能量,這是一種純粹的量子效應,交換能量通常在
化學鍵結中扮演關鍵角色, r12 代表二個電子間的距離。。
Hartree-Fock 方法最複雜的地方在於怎麼去最佳化Eq. (4) 中的能量期望值,或換
句話說,怎麼去最佳化 Eq.(3)-(7) 中的分子軌域 (MO)。經過複雜的理論推導,我們可
以證明最佳化的分子軌域要滿足以下的 Hartree-Fock Equation
Fˆ(1)φi (1) = ε iφi (1) (8)
其中F 稱為 Fock operator,它是一個很複雜的單電子運算元,其中包含了:(a) 計算在
φi 軌域電子的動能以及它與所有原子核間的作用力的運算元 (b) 計算在 φi 軌域電子
所參與的電子間庫倫作用力的運算元 (Coulomb operator) (c) 計算在 φi 軌域電子所參
與的電子間交換作用力的運算元 (exchange operator),ε
i 為分子軌域能量。 雖然 Eq. (8)
4
看起來像一個單純的 eigenfunction-eigenvalue 問題,但是由於Fock operator 中的庫倫
及交換運算元本身定義上需要先知道 eigenfunctions φ
i,因此實際上在解 Eq. (8) 時需
要先給定一組近似的起始分子軌域 (initial guess),然後以自相吻合 (self-consistent field,
SCF) 的反覆計算來求解,直到 Hartree-Fock (HF) 能量收斂到沒有顯著變化為止。值
得注意的是在HF 理論中軌域能量的總和並不等於電子的總能量,因為在分子軌域能量
中庫倫及交換作用力會被重複計算,我們可證明 Hartree-Fock 方法所得到的總能量也
可寫成:
( ) NN
n
i
n
j
ij ij
n
i
E = Σ i −Σ Σ J −K +V
=1 =1 =1
HF 2 ε 2 (9)
以上的運算都是針對分子軌域;然而,我們應該用什麼數學方式來表示分子軌域呢?最
常用的方法是將每一個MO 可寫成分子中所有原子軌域 (AO) 的線性組合,假設共使
用K 個AO,則
Σ=
=
K
j
i c ji f j
1
φ (10)
由此,對MO 最佳化的抽象概念轉變成了找出最洽當的展開係數 (cji) 的具體目標。這
K 個AO 通常稱為計算使用之基底函數 (basis functions or basis set)。Roothaan 在1951
年提出以此 LCAO-MO 為架構求解 Hartree-Fock Equation 的方法,基本上就是將此複
雜的eigenvalue-eigenfunction 的問題以線性代數的方法來有效的處理,Eq. (8) 可以用
矩陣方法很簡單的表示成
FC = SCε (11)
F, S 分別稱為 Fock matrix 及 overlap matrix,定義為
( ) ( )⎥⎦⎤⎢⎣= + ⎡−⎥⎥⎦⎤⎢⎢⎣⎡−+ ⎟⎟⎠⎞⎜⎜⎝⎛= −∇ −=
Σ Σ
∫ ∫
∫ Σ Σ Σ Σ
∫
= =
= = =
H P rs tu ru ts
d d
r
d d f f f f
r
f f f f
f d e c c
r
Z e
m
f
F f F f d
K
t
K
u
rs tu
r s t u r u t s
uj
K
t
K
u
n
j
s tj
A A
A
e
r
rs r s
|
2
| 1
2 (1) (1) (2) (2) (1) (1) (2) (2)
4πε(1)
4πε2
(1) *
ˆ1 1
core
1 2
12
* *
1 2
12
* *
*
1 1
/ 2
1
*
0
2
1
0 1,
2
2
1
2
τ τ τ τ
τ
τ
h
(12)
Srs = ∫ fr fs dτ (13)
5
而 C 稱為 coefficient matrix, 也就是由 (10) 式中的 cji 所組成的,F, C, S 都是 K × K
的矩陣。ε 是一個K × K 的對角線矩陣,對角線上的值就是 (8) 式中的分子軌域能量,
(rs | tu) 及 (ru | ts) 為Coulomb 及 exchange 雙電子積分的常用縮寫。 Eqs. (11-12) 的
推導與詳細的解法超出本章所要介紹的範圍故在此省略。Eq.(12) 中的Ptu 稱為density
matrix elements,因為我們可證明在 HF 理論中分子內的電荷分布與P 有密切的關係:
( ) Σ Σ
= =
=
K
r
K
s
r Prs fr fs
1 1
ρ * (14)
而HF 能量也可以寫成
( ) NN
K
r
K
s
E Prs Frs Hrs V
2
1
1 1
core
HF = Σ Σ + +
= =
(15)
由 (12) 式我們可以看出每一個 Fock matrix element Frs 的計算量大約與 K 平方成正
比,而Fock matrix 總共有K 平方個Frs,因此HF 方法的計算量大約跟K 的四次方成
正比,所以HF 方法的計算量會隨著化學系統或基底函數的增大而很快的上升。然而,
因為 HF 方法只考慮了電子間的平均作用力,通常無法得到非常準確的分子能量,但
HF 方法在許多情況下仍然能夠提供一些很有用的定性預測與最佳化的分子軌域,對於
穩定的分子,HF 方法通常也能預測出非常準確的分子結構。若要更進一步得到更準確
的能量,我們需要利用更複雜的理論來考慮到電子間瞬間的作用力,或說是要計算所謂
的電子相關能量 (electron correlation energy)。
(2) Configuration Interaction (CI)
在 HF 計算中Eq. (11) 的解讓我們得到 K 個 MO,在 (3) 中使用最低能量的 n
個 MO 可以讓我們得到所謂的HF wavefunction。也就是說在 HF 方法中我們將2n 個
電子指定於最低能量的n 個軌域,這種將電子指定於軌域的概念稱為一種組態或
configuration。一般而言,我們可以利用這K 個MO 產生許多種不同的configurations
或 slater determinants,而HF 理論只使用了一個最低能量的 configuration。因此,一種
最直接的改進HF 理論的方法就是將波函數寫成所有configuration 所對應之 slater
determinants 的線性組合:
Ψ = Ψ +Σ Ψ + Σ Ψ + Σ Ψ + Σ Ψpqrs +L
ijkl
p q r s
i j k l
pqrs
ijkl
pqr
ijk
p q r
i j k
pqr
ijk
pq
ij
p q
i j
pq
ij
i p
p
i
p
c ci c c c
, , ,
, , ,
, ,
, ,
,
, ,
CI 0 HF (16)
其中 i, j 等代表在 HF 理論中填滿的軌域 (occupied orbitals),而 p, q 等代表HF 理論
中的空軌域 (virtual orbitals)。因此 p
Ψi 代表將一個原來在第 i 個填滿軌域的電子放到
6
第 p 個空軌域所形成的所謂singly excited slater determinants,而 p q
i j
,
Ψ, 代表將原來在第
i 及 j 填滿軌域的二個電子放到第 p 及 q 二個空軌域所形成的所謂doubly excited slater
determinants,依此類推。每一個行列式前的係數則由 variational principle 將能量最小
化來決定。若 (16) 式中包含所有可能的組態則稱為 Full Configuration Interaction (FCI)
理論,我們可以證明 FCI 理論相當於求解完整的薛丁格方程式。實際上由於 FCI 理論
所產生的configuration 數量過於龐大,對於一般的分子而言遠遠超過電腦的負荷量,
因此在 (16)中必須做一些省略 (limited-CI or truncated CI),最常用的 CI-SD 方法就是
只使用 (16) 中的前三項而忽略其他 highly excited configurations,更精確但極為費時的
CI-SDTQ 則是使用(16) 中的前五項。單純的CI-SD 方法現在很少使用在化學反應系統
上,因為通常truncated CI 所得到的相對能量並不可靠,一種較新的 quadratic CI 方法,
如QCISD 或 QCISD(T) 方法則可以大幅改善相對能量的問題。
(3) M鴏ler-Plesset Purterbation Theory (MPn)
早在1934 年M鴏ler 及Plesset 提出了一種計算分子系統的方法,這種方法是根據
perturbation theory,就是將系統的 Hamiltonian 分成二部份,一部分是可直接求解的
H°,另一部份是難以求解的 perturbation 項 H΄Hˆ= Hˆ° + Hˆ' (17)
在M鴏ler-Plesset 的方法中,H° 為單電子Fock operator 之和:
Σ=
° =
n
i
H f i
2
1
ˆˆ( ) (18)
因此,由 (3), (8)
HF
1
HF ε2 ˆΨ ⎟⎟⎠⎞⎜⎜⎝⎛= Ψ ° Σ=
n
i
H i (19)
因此,在 unperturbed 系統中,能量 (zeroth-order energy) 為 MO 能量之和。此方法中
第一階基態能量校正量為
E ψ*Hˆ' ψdτ HF Hˆ' HF dτ
) 0 (
0
) 0 (
0
) 1(0 = ∫ = ∫Ψ Ψ (20)
所以經過第一階校正後的總能量可以寫成
7
( ) HF ( ) HF HF
) 0 (
0
) 0 (
0
) 1(0
) 0 (
0 E + E = ∫ψ* Hˆ° + Hˆ' ψdτ = ∫Ψ Hˆ° + Hˆ' Ψ dτ = E (21)
因此,經過第一階校正我們又回到 HF 理論的能量,所以額外的校正是要由第二階開
始。由perturbation theory 及所謂的 Condon-Slater rules 我們可以得到:
( ) ( ) Σ ΣΣΣΣ
∫
+ −+
−=
−=
m≠ m a b i j i j a b
m ij ab ia jb
E E
ψH ψd
E
εεεε* ˆ' | | 2
0
) 0 ( ) 0 (
0
2 ) 0 (
0
(0)
) 2 (
0
τ
(22)
其中i, j 代表不同的occupied spin orbitals,a, b 代表不同的virtual spin orbitals。上式中
的) 0 (
mψ為有使用virtual orbitals 的較高能量之slater determinants,在第二階能量校正中,
只有使用二個virtual orbitals 的 ) 0 (
mψ(doubly excited determinants) 有貢獻。這種加入第
二階校正的方法稱為 MP2:
) 2 (
0 EMP2 = EHF + E (23)
MP2 的計算量大約與系統大小的五次方成正比。MP2 是一種兼顧準確性與計算效率的
一種非常重要的方法,在第二階校正中此方法將大部分的電子相關效應考慮進去,使得
MP2 方法可以遠較HF 方法更準確的計算出大部分化學系統的相對能量及更可靠的分
子結構,雖然MP2 的計算量通常明顯的較 HF 為高,但現在一般計算化學用的電腦可
以使用MP2 理論很容易的處理大約二十個原子以下的系統。依照pertubation theory 我
們可以進一步加入第三、第四階的校正,而所得到的方法稱為 MP3 及 MP4。一般而
言,MP3 並不能有系統的增加 MP2 的準確度, MP4 方法通常可以再顯著提高準確
度,但其計算量會顯著的增加,而且大約與系統大小的六或七次方成正比,不容易應用
於較大的系統中。
(3) Coupled-Cluster Method (CC)
Coupled cluster 是早在1960 及1970 年代間發展出來的準確量子化學方法,由於它
的計算量很大,一直到了最近十年來才被廣泛的使用在一般的化學系統的計算。CC 方
法的基本方程式為
HF
ˆΨ = eTΦ (24)
8
= + + + +L
3 !
ˆ2 !
ˆˆ1
ˆT 2 T 3 eT T (25)
TˆTˆTˆTˆTˆn = 1 + 2 + 3 +L+ (26)
其中Tˆn是由HF波函數產生所有n-electron excitation之slater determinants的運算元,例
如:
Σ
Σ
Ψ = Ψ
Ψ = Ψ
a b
i j
ab
ij
ab
ij
i a
a
i
a
i
T t
T t
,
,
2 HF
,
1 HF
ˆˆ(27)
方程式 (24) 可被證明是 FCI 的另一種寫法,而各項係數 t 則須由滿足 (24) 式的前提
下經由複雜的推演而導出的一系列龐大的非線性方程組中求反覆求解。如同 CI 理論,
完整的 CC 理論對於一般化學分子而言計算量過於龐大,因此必須做一些簡化。通常
的簡化方法是只保留 (26)式中的頭幾項而忽略高次激發項。比如說,若只用 T2 則稱為
CCD 理論,若保留了 T1 及T2,則稱為 CCSD 理論。其實T3 項對於達到高準確度而言
非常重要,只是若同時保留T1,T2,及T3 則計算量通常會變的難以負荷,因此一種折
衷的辦法是不直接使用T3 而是利用 perturbation theory 等方法去估計T3 的貢獻,如現
在常被用來當作高準確度量化計算之標準方法的 CCSD(T) 理論就是屬於這種類型。一
般認為 QCISD(T) 及 CCSD(T) 方法是當今常用的量化計算理論中最可靠的二種方
法,不過它們必須搭配適當的基底含數才能完全發揮功能。
(4) 基底函數 (Basis Set)
在上述的量化計算方法中都使用到分子軌域 (MO) 的觀念,式 (10) 中所用到用來展開
分子軌域的原子軌域 (AO) 稱為量化計算之基底函數 (basis functions),傳統上用來做
多電子原子及雙原子分子計算的基底函數是所謂的 Slater-type orbital (STO):
1 / 0 ( , )
STO θ φ m
l
f = N rn−e−Zr a Y (28)
STO 是類氫原子軌域的簡化,此處的 Z 稱為 orbital exponent,N 為 normalization
constant,STO 可以正確的描述在原子核附近波函數的行為。在 HF 計算中每一個 MO
是以數個 STO 線性組合而成。但在多原子分子的計算時STO 非常沒有效率,因此Boys
等人在 1950 年代提出使用所謂的 Gaussian-type functions (GTO):
9
2
GTO
f = N xi y j zk e−Zr (29)
其中i, j, k 為零或正整數, Z > 0 稱為 orbital exponent,N 為 GTO 之normalization
constant。當i + j + k = 0 時稱為 s-type GTO,當i + j + k = 1 時稱為 p-type GTO,當i +
j + k = 2 時稱為 d-type GTO,依此類推。由 (29) 在同一個Z 值下我們可以得到六種
d-type GTO,通常的作法是將其線性組合成類似實數3d AO (dxy, dxz, dyz, dx2-y2, dz2) 的
五個 GTO 而省略具有s 對稱性的一個 GTO。 式 (29) 之 GTO 也稱為所謂的
Cartesian Gaussian,其中AO在角度上的變化是以簡單的x, y, z 函數來取代複雜spherical
harmonic 函數,在exponential 項上是使用 r 平方而非 STO 中的 r 一次方。上述之
STO 及GTO 都是以原子核為中心之AO。使用 GTO 可以大幅簡化雙電子積分的計
算,因為以二個不同原子為中心的 GTO 的乘積等於另一個以這二個原子之間的點為中
心的 GTO。但為了能如 STO 一樣正確描述原子核附近波函數的行為,通常我們需要
將數個 GTO 做線性組合成一個行為上類似 STO 函數的 contracted Gaussian-type
function (CGTF):
=Σ
l
fCGTF dl gl (30)
其中gl 為以同一個原子為中心的數個 Cartesian Gaussian (29),但具有不同的 exponents
(Z),dl 為展開係數,gl 也常被稱為所謂的 primitive Gaussians。用來描述一種原子的
所有 CGTF 稱為基底函數或 basis set。 對每一個原子而言,如果 CGTF 的數目與其
在週期表中同週期原子可用之原子軌域數相同,則稱為 minimal basis set。比如說對碳
原子而言,minimal basis set 會包含一個 s-type 的 CGTF 描述1s orbital,另一個 s-type
的 CGTF 描述2s orbital,另一組(3 個) p-type 的CGTF 描述2p orbitals。量化歷史發展
上非常有名的STO-3G basis set 就是屬於這種minimal basis set,其中每一個CGTF 是
利用三個GTO 線性組合而成來模擬一個 STO AO。
Minimal basis set 所得的計算結果通常最好也只能算是做到定性的預測,要進一步
增加準確度需要增加基底函數的量。所謂 double-zeta (DZ) basis set 是指對每一個可用
的原子軌域而言使用二個CGTF 來描述,使得計算上用到的基底函數之數目加倍,比如
說Dunning 及 Huzinaga 的 D95 basis set 就是屬於此類型。DZ basis set 會使計算量大
量增加,一種折衷的辦法是只將價軌域改成 DZ,內層軌域維持minimal basis set,因為
內層軌域的貢獻通常在計算相對能量時會抵消掉。此種基底函數稱為 split-valence
double-zeta basis set,如常見的 D95V, 3-21G, 6-31G 等。
常見的3-21G, 6-31G 等基底函數稱為 Pople-type basis set。在 3-21G 中,每一個
內層電子 (core electron) orbital 是由三個 primitive Gaussians 所組成的一個 CGTF 來代
表,每一個價電子軌域則是由二個CGTF 來代表,其中一個CGTF 是由二個primitive
10
Gaussians 所組成的,而另一個CGTF 則是一個 exponent 絕對值最小的 uncontracted
GTO。在 6-31G basis set 中情況也類似,每一個內層電子 (core electron) orbital 是由六
個 primitive Gaussians 所組成的一個 CGTF 來代表,每一個價電子軌域則是由二個
CGTF 來代表,其中一個CGTF 是由三個primitive Gaussians 所組成的,而另一個CGTF
則是一個 exponent 絕對值最小的 uncontracted GTO。
現在使用的基底函數同常會加上所謂的 polarization functions,也就是具有比價軌
域更高的角動量量子數之 AO,比如說在所謂的 6-31G* 或 6-31G(d) 的基底函數中,
對於所有第二週期 (Li-Ne) 及第三週期 (Na-Ar) 的原子都加上了一組 uncontracted
d-type GTO。加入polarization functions 的目的是在分子的計算中較容易將電子的密度
朝鍵結的方向極化,得到比較可靠的結構與能量。在所謂的 6-31G** 或 6-31G(d,p) 的
基底函數中,對於氫及氦原子也加入了一組 p-type polarization GTO。在更精確的計算
中,我們需要用到更大的基底函數,比如說 6-311G basis set 是一個 valence triple-zeta
的基底函數,對於所有第二週期原子每一個價電子軌域則是由三個CGTF 來代表,其
中一個CGTF 是由三個primitive Gaussians 所組成的,而另二個CGTF 則是由二個
exponent 較小的 uncontracted GTO 所組成。同樣的,6-311G 也可加入極化函數形成如
6-311G** 或 6-311G(d,p) 等基底函數。有時候多加入幾組極化函數對於能量及一些性
質的預測會更為準確,例如 6-311G(2df,2pd) basis set 代表對第二週期及以後的原子加
入二組d-type 及一組 f-type 極化函數,並且對第一週期原子加入二組p-type 及一組
d-type 極化函數。在研究陰離子、凡德瓦爾作用力、反應過渡態時由於電子雲分佈的範
圍比較廣,我們常需要使用到涵蓋空間較大的分子軌域,因此我們需要加入一些所謂的
diffuse functions,也就是 orbital exponents 的絕對值比較小的基底函數,例如
6-31+G*、6-311+G* 代表在 6-31G* 或 6-311G* basis set 中再加入一組 s- 及 一組
p-type 的diffuse functions,而 6-31++G* 或 6-311++G* 則代表對氫及氦也加入一組
s-type 的diffuse functions。通常對第一週期原子加入 diffuse functions 的效用並不明顯。
Dunning 等人自從 1989 年起發展了另外一個系列的基底函數,稱為
correlation-consistent (cc) basis set (cc-pVnZ, n = D, T, Q, 5, 6),他們著重的地方在於高階
correlation energy 的計算,以及外插至 complete basis set (CBS) limit 的收斂情形。在這
些 basis sets 中, polarization function (p) 是內含的,VDZ 代表 valence double-zeta,
VTZ 代表 valence triple-zeta,依此類推。對第二週期及之後的原子而言DZ 中含有d
polarization functions, TZ 中含有d, f polarization functions,QZ 中含有d, f, g
polarization functions,依此類推。對第一週期的原子而言DZ 中含有p polarization
functions, TZ 中含有p, d polarization functions,QZ 中含有p, d, f polarization
functions,依此類推。這一類基底函數的 diffuse functions 可藉由加上 aug- (augmented)
的字頭來指定,比如說 aug-cc-pVDZ 是指原來的cc-pVDZ basis set 中再在加上一組 s,
p, d diffuse functions.。在使用高階理論 (如MP4、QCISD(T)、CCSD(T) 等) 計算準確
相對能量或要外插到 CBS limit 時,correlation-consistent basis sets 是較好的選擇。
11
通常進行量化計算時需要同時指定理論方法以及基底函數,符號上一般是以
theory/basis 的方式表示,例如 HF/3-21G,MP2/6-31+G**,CCSD(T)/aug-cc-pVTZ 等。
以計算相對能量而言 HF 方法所得的結果對基底函數的大小不太敏感,使用valence
double-zeta 以上的 basis set 通常沒有什麼幫助;然而計算 correlation energy 時,basis
set 的品質就非常重要,在 MP2 的計算中,由double-zeta 到 triple-zeta 通常相對能量
會有顯著改進,而 polarization functions 也是得到準確能量所必需的。對於更高階的理
論如 CCSD(T), QCISD(T) 等,一定要搭配上很好的基底函數 (如 aug-cc-pVTZ 等) 才
可以充分發揮效能。
(5) Density Functional Theory
以薛丁格方程式求波函數是一項艱鉅的工作,在n 個電子的系統中波函數包含了
3n 個空間座標以及 n 個自旋座標,並且必須是 antisymmetric ,而這個極為複雜的波
函數似乎也沒有很直接的物理意義,只能靠量子力學中的 operator 去擷取所需的資
訊。換句話說,複雜的電子波函數似乎包含了太多一般化學家所不需要的資訊。 那麼
有沒有其他較簡單的函數足以描述量子化學系統的性質且較容易計算出來呢?一個化
學系統的電子密度分佈似乎滿足這個要求,因為在 Born-Oppenheimer 假設下電子密度
決定了整個化學系統所有的電荷分佈,而通常電荷分佈決定了化學分子的許多重要的性
質。的確,在 1964 年Pierre Hohenberg 以及 Walter Kohn 證明了分子基態的能量及其
他所有性質可以由基態的電子密度分佈唯一決定,或者說基態的性質是電子密度分佈的
泛涵 (functional),此處的functional 是指函數的函數,也就是說分子性質由整個密度分
佈函數來決定。 因此理論上電子密度函數的確在很多情況下可以取代波函數來描述分
子的各項性質,但問題是 Hohenberg-Kohn 的理論並沒有告訴我們怎麼樣不使用波函數
去求電子密度,也沒有告訴我們怎麼由電子密度去求分子的能量。
1965 年時 Kohn 及 Sham 提出一個實際上使用Hohenberg-Kohn 的理論的方法,
一般稱為所謂的 Kohn-Sham (KS) method。這個方法是假設一個虛擬的系統,電子的數
目與實際的系統一樣,但電子之間並沒有作用力,但每個電子都感受到一種位能,或稱
為 external potential。本來電子感受到的位能是來自於電子之間以及電子與原子核之間
的作用力,但在此虛擬系統裡,電子所感受到的external potential 是一種很特別的型態,
它的定義是在這種 external potential 下,虛擬系統的電子密度函數會等於真實系統下的
電子密度函數。KS method 的功能是提供了我們在實際操作上得到近似的電子密度以及
能量的一種方便方法。因為在虛擬系統中電子間沒有作用力,電子的波函數是一個
slater determinant,其中的 orbital 滿足
ˆKS KS KS KS hi θi =ε i θi (31)
12
hi i us (ri )
2
ˆKS = −1 ∇2 + (32)
其中 us 是上述所謂的 external potential,θKS 稱為 density orbital 或 Kohn-Sham (KS)
orbital。由此,我們可以證明系統的電子密度可以表示成
Σ=
=
n
i
i
1
KS 2 ρ θ (33)
而此虛擬系統中電子的動能可以寫成
Σ=
= −∇
n
i
Ks i i
1
2 KS
1
KS(1) | | (1)
2
1 θ θ (34)
電子與原子核間的位能可寫成
1
1
(r1) dr
r
UNe = −ΣZ ∫
α α
α
ρ
(35)
其中 α 代表不同的原子核,Zα 代表核電荷,r1α 代表空間中任一點到原子核 α 的距
離,而電子之間的古典靜電位能可以寫成
1 2
12
(r1) (r2 ) r r
2
1 d d
r
Uee = ∫ ∫
ρ ρ
(36)
其中的1/2 是避免重複計算電荷間的作用力。 以上這些能量的計算都很直接,只要有
了上述的 Kohn-Sham orbitals 以及密度分佈函數 ρ(r) 就可以很容易的得到。然而,式
(34) 與 (36) 有一些其他的問題。(34)是虛擬系統中的動能,真實系統的動能必然有些
不同;(36) 是古典靜電能量,並不包含量子力學中的電子間的 exchange 及 correlation
energy。我們若將這些校正量合稱為 exchange-correlation energy functional, Exc,則系統
的總能量可以寫成
E =UNe + Ks +Uee + Exc (37)
Hohenberg and Kohn 也證明了如果我們知道真正的density functional,則真正的電子密
度會最小化系統的能量,這稱為 Hohenberg-Kohn variational theorem。因此,我們可以
藉著改變電子密度分佈函數來最小化式 (37) 所得到的能量,以求得系統的電子密度與
能量。由式 (33) 我們知道電子密度可由 KS orbitals 求出,因此我們可以藉由改變 KS
orbitals 來改變電子密度分佈函數。最佳化 (37) 所得到的能量我們可以得到一個類似
Hartree-Fock Equation 的方程式,稱為Kohn-Sham Equation
13
ˆKS(1) KS(1) KS(1) KS(1)
h θi = ε i θi (38)
其中 KS Hamiltonian 可表示成
⎥⎥⎦⎤⎢⎢⎣⎡= −∇ −Σ + ∫ +
α α
α ρ
(r ) r (1)
2
ˆ(1) 1 2
12
2
1
2
1
KS
d vxc
r r
h Z (39)
其中 vxc 是Exc 對密度分佈函數的微分
δ(r)
(r) δ[ (r)]
ρ
xc ρ
xc
v = E (40)
式 (38) 可用類似 HF 理論中的 SCF 方法求解。如同 molecular orbital,KS orbital 可
以寫成是atomic orbital 或基底函數的線性組合,因此求解 (38) 也等於求構成KS
orbital 最佳的基底函數線性組合。到目前為止 DFT 理論並沒有包含任何的近似,所以
理論上DFT 可以得到完全正確的結果,但不幸的是我們並不知道 Exc 真正的數學型態,
只能用不同的方法去近似它,不同的 DFT 方法的差異就在於使用不同的 exchangecorrelation
energy functional。如果Exc 只跟密度函數本身有關,我們稱為 Local (Spin)
Density Approximation (LDA or LSDA),現在化學家常使用的DFT 方法中其Exc 也跟密
度函數的 gradient 有關,也稱為 generalized-gradient approximation (GGA)。 一般通常
將 Exc 分開成 exchange 與 correlation 項
Exc = Ex + Ec (41)
而DFT 方法通常是由所使用的 exchange 及 correlation functionals 來命名,比如說
BLYP 方法是指使用 Becke 在 1988 年發表的 exchange functional (B 或 B88) 以及
使用 Lee-Yang-Parr 三個人在1988 年所發表的 correlation functional,而 BPW91 方法
中則是使用上述的exchange functional 以及 Perdew-Wang 二人在 1986 年發表的
correlation functional。近年來的研究發現,如果在 exchange functional 中混入一定比例
的類似Hartree-Fock 的exchange energy 但使用的是KS orbitals
Σ Σ
= =
= −n
i
n
j
x i j r j i
E
1 1
KS KS
12
HF KS(1) KS(2) | 1 | (1) (2)
4
1 θ θ θ θ (42)
則在化學計算上的準確度可以大幅提升,此種方法稱為 Hybrid DFT,其中最有名的莫
過於在1993-1994 年所發展出的 B3LYP 方法
14
LSDA LYP LSDA
0
HF B88
0
B3LYP (1 ) (1 ) Exc = a Ex + axEx + −a −ax Ex + acEc + −ac Ec (43)
由與實驗值對照所得之最佳化係數為 a0 = 0.20, ax = 0.72, ac = 0.81。B3LYP 方法是近十
餘年來最受歡迎的 DFT 方法,因為幾乎在各種型態分子中的不同性質都可以由此方法
得到可靠的結果,許多人甚至將DFT 與 B3LYP 劃上等號。B3LYP 最常搭配的基底函
數包括6-31+G(d,p), 6-311+G(2d,p), 6-311+G(3df,2p) 等。然而,最近幾年來發展出了一
些新的 hybrid functionals 包括 mPW1PW91, B1B95, B98, BMK, M05, M06, M08 等其
準確度在很多方面已經可以超過 B3LYP。Hybrid DFT 的計算量大約與系統大小的四次
方成正比,與 Hartree-Fock 方法類似,但現代的hybrid DFT 的準確度可以與MP2 方
法相當甚至更好,一般認為DFT 在未來的發展上應該還有很大的改進的空間,最近的
研究也顯示,DFT 可以與wavefunction 的方法相結合而得到很準確又有效率的理論方
法。
胡維平
國立中正大學
化學暨生物化學系
?Copyright 2009
No comments:
Post a Comment