1
量化計算方法簡介
很多人對量化計算很排斥，因為不少學化學的人覺得量化實在太難了。量子
化學的確不容易學的精通，但就如做實驗一樣，你不用完全了解實驗背後的複雜
化學理論與儀器原理，只要你能掌握好實驗條件與基本的操作技巧，你一樣可以
得到有意義的數據。在現代的量化計算中，由於真正繁瑣複雜的工作都是由電腦
軟硬體來完成，你其實也只要對電腦的基本操作熟悉，對量子化學理論以及分子
結構有基本的認識，就可以從事有意義的計算工作。 當然，對量子化學有比較
清楚的認識，對於計算方法的了解以及對計算結果的解讀都會有很大的幫助。以
下我針對基礎的量化計算所需具備的量化理論做一個簡要的介紹。
1. 薛丁格方程式 (Schrödinger Equation)
對於大部分的基態 (electronic ground-state) 獨立的 (free) 化學分子，理論上
我們可以藉著解下列與時間無關的薛丁格方程式 (time-independent Schrödinger
Equation) 而得到分子的所有性質：
0 0 0 Hˆψ = E ψ (1)
其中Hˆ 為分子系統的總能量運算元 operator，其中包括電子的動能 (Ke)，電子
間的位能 (Vee)，電子與原子核間的位能 (VNe)，以及原子核間的位能 (VNN)。
Hˆ = Kˆ e +Vˆee +VˆNe +VˆNN (2)
而 ψ0 為基態波函數，它是一個包含了所有電子位置座標以及自旋座標的複雜函
數，E0 則為基態能量。方程式 (2) 隱含著所謂的 Born-Oppenheimer 假設，也
就是說原子核的運動可以和電子的運動分開處理，因次在 (1) 式中我們先假設原
子核是靜止的。 此外，量子力學理論也要求波函數必須是 antisymmetric，也就
是說任意調換二個電子的所有座標波函數必須變號，而且理論上波函數也應該是
total spin operator (所有電子自旋角動量的向量和) 的 eigenfunction:
ˆ ( 1) , 0, 1/2, 1, 3/2, ... 0
2
0
S 2ψ = S S + h ψ S = (3)
2
對於多電子的化學系統，求解 (1) 式是非常困難的工作，因此在實際的量化計
算中，我們通常以各種近似的理論利用數值方法配合電腦的大量高速運算來求得
(1) 式的近似解。
2. Hartree-Fock (HF) Method
HF 方法是最基本的一種量化計算理論，大部分更精確的方法都是 HF 方法
的延伸。HF 方法概念上並不困難，主要就是以行列式來近似分子的波函數，以
closed-shell 的分子為例：
(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
HF
n n n n n n
n
n n
n n
n n
φ α φ β φ α φ β φ α φ β
φ α φ β φ α φ β φ α φ β
φ α φ β φ α φ β φ α φ β
LL
M
M
LLL
LLL
Ψ = (4)
其中φn 為被電子填滿的分子軌域`(MO), 1...2n 代表不同的電子，此種行列式也被
稱為 slater determinants。這其實跟量化課程中所學到分子近似波函數沒什麼不
同，以 H2 分子為例，電子的近似波函數常寫成：
(2) (2)
(1) (1)
2
1
1 1s
1 1s
H2 φ α φ β
φ α φ β
σ σ
σ σ
s
Ψ = s (5)
這其實就是一種 HF 波函數，但在真正的計算化學中，分子軌域不再侷限於簡
單的類氫原子價軌域 (1s, 2s, 2p, ... etc.)的線性組合，而通常是由許多個具有相同
對稱性的基底函數 (basis functions) 組合而成。此外，所有基底函數線性組合的
係數是利用 variational theory 將能量的期望值最小化而得到：
H dτ
cri
Min ∫ΨHF ˆΨHF (6)
其中 cri 是指第 i 個 MO 中第 r 個基底函數的展開係數。此能量最小化的過程
非常複雜，最主要的計算時間也都花在這個步驟，簡單來說，(6) 式無法直接求
解，因為過程中需要分子軌域的函數，但在 (6) 式解出前我們並不知道分子軌
域的線性組合資料。因此 (6) 式的求解過程我們通常稱為 self-consistent field
(SCF)方法，首先電腦程式會先產生一組近似的分子軌域，稱為 initial guess，然
後利用 SCF 的方法逐漸收斂到產生 (6 )式中最低能量的一組 MO，此時的能量
稱為 HF energy。在 HF 計算中一共會得到的MO 數目與所使用的基底函數的數
3
目 (K) 相同，假設系統中有2n 個電子，其中最低能量的n 個MO 為occupied
MOs，其他則為unoccupied MOs 或稱為 virtual orbitals，通常virtual orbitals 的
數目會遠大於 occupied MO 的數目。 由 (2) 式可以看出HF 能量包括電子的動
能 (正值) 及位能 (通常為負值) 以及原子核間的位能 (正值)，而能量的零點是
所有電子與原子核分離到無窮遠而且處於靜止狀態。 通常量化計算所得的能量
是以 hartree (h) 來表示：
1 hartree = 627.5095 kcal/mol (7)
依此單位氫原子在基態時的能量為 −0.500 h.
最後有幾件是值得在此說明，首先 (1), (2) 式是根據所謂 non-relativistic 的
量化理論，也就是說沒有考慮相對論效應。西元 1928 年 Dirac 發現自旋其實
是一種相對論效應，但完整考慮相對論效應的量子力學理論在多電子系統下非常
難以計算，因此電子自旋的效應在一般 non-relativistic 的量化理論中是以外加的
方式考慮，幸好這對於前三週期的化學元素並不會產生顯著的誤差，但此理論並
無法處理我們在量化課程中所學到的 spin-orbital coupling 等直接與自旋相關的
效應，以及高速運動時電子質量增加的現象。此外， HF 理論只是一種求解 (1)
式的近似方法，此方法雖然可以滿足 (對 closed-shell 分子而言) 波函數
antisymmetric 的要求，但此理論只有在單電子系統中是完全正確的，這是因為
我們將波函數侷限在單一行列式的型態中，有效而言這只考慮了電子間的平均作
用力，但忽略了電子間的瞬間作用力。HF 方法未考慮到的效應所產生的能量貢
獻一般統稱為 correlation energy (CE)。 在許多化學性質的計算中，correlation
energy 都會有很重要的貢獻，因此現代利用波函數的量化計算通常是以 HF 方
法為基礎再加上不同程度的近似去估計 correlation energy 的效應，如接下來要
提到的 MP2，CCSD 等方法。
3. Møller-Plesset Purterbation Theory (MPn)
早在1934 年M鴏ler 及Plesset 提出了一種計算分子系統的方法，這種方法
是根據 perturbation theory，就是將系統的 Hamiltonian 分成二部份，一部分是
可直接求解的unperturbed H°，
(0)
0
(0)
0
(0)
Hˆ° ψ0 = E ψ(8)
其中(0)
ψ0 和 (0)
E0 分別稱為是zeroth-order 基態波函數及基態能量。
另一部份是難以求解的 perturbation 項 H΄，因此
Hˆ= Hˆ° + Hˆ' (9)
4
我們希望將基態能量寫成zeroth-order 的基態能量加上各階校正量
= + + + + (4) +L
0
(3)
0
(2)
0
(1)
0
(0)
E0 E0 E E E E (10)
簡單來說Perturbation theory 的精神在於利用zeroth-order 的所有波函數及能量
以及H΄推導出將各階能量的校正量。在M鴏ler-Plesset 的方法中，H° 為單電子
Fock operator 之和，詳細的理論推導在此先省略。
我們可以證明在M鴏ler-Plesset 的方法中經過第一階校正後的總能量可以寫成
( ) ( ) HF HF HF
) 0 (
0
) 0 (
0
) 1(0
) 0 (
0 E + E = ∫ψ* Hˆ° + Hˆ' ψdτ = ∫Ψ Hˆ° + Hˆ' Ψ dτ = E (11)
因此，經過第一階校正我們又回到 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
τ
(12)
其中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 (13)
MP2 的計算量大約與系統大小的五次方成正比。MP2 是一種兼顧準確性與
計算效率的一種非常重要的方法，在第二階校正中此方法將大部分的電子相關效
應考慮進去，使得MP2 方法可以遠較HF 方法更準確的計算出大部分化學系統
的相對能量及更可靠的分子結構。雖然MP2 的計算量通常明顯的較 HF 為高，
但現在一般計算化學用的電腦可以使用MP2 理論很容易的處理大約二十個原子
以下的系統。依照pertubation theory 我們可以進一步加入第三、第四階的校正，
而所得到的方法稱為 MP3 及 MP4。一般而言，MP3 並不能有系統的增加 MP2
的準確度， MP4 方法通常可以再顯著提高準確度，但其計算量會顯著的增加，
5
而且大約與系統大小的六 (MP4SDQ) 或七次方 (MP4SDTQ) 成正比，不容易應
用於較大的系統中。
4. Configuration Interactions (CI)
在 HF 計算中讓我們得到 K 個 MO，在 (4) 式中使用最低能量的 n 個
MO 可以讓我們得到所謂的HF wavefunction。也就是說在 HF 方法中我們將2n
個電子指定於最低能量的n 個軌域，這種將電子指定於軌域的概念稱為一種組態
或configuration。一般而言，我們可以利用這K 個MO 產生許多種不同的
configurations 或 slater determinants，而HF 理論只使用了一個最低能量的
configuration。因此，一種最直接的改進HF 理論的方法就是將波函數寫成所有
configurations 所對應之 slater determinants 的線性組合：
Ψ + Ψ +L
Ψ + Ψ +
Ψ = Ψ +
Σ Σ
Σ Σ
pqrs
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
i
c c
c c
c
, , ,
, , ,
, ,
, ,
,
, ,
CI 0 HF
(14)
其中 i, j 等代表在 HF 理論中填滿的軌域 (occupied orbitals)，而 p, q 等代表
HF 理論中的空軌域 (virtual orbitals)。因此 p
Ψi 代表將一個原來在第 i 個填滿
軌域的電子放到第 p 個空軌域所形成的所謂singly excited slater determinants，
而 p q
i j
,
Ψ, 代表將原來在第 i及 j填滿軌域的二個電子放到第 p及 q 二個空軌域
所形成的所謂doubly excited slater determinants，依此類推。每一個行列式前的係
數則由 variational principle 將能量最小化來決定。若 (14) 式中包含所有可能的
組態則稱為 Full Configuration Interaction (FCI) 理論，我們可以證明 FCI 理論相
當於求解完整的薛丁格方程式。實際上由於 FCI 理論所產生的configuration 數
量過於龐大，對於一般的分子而言遠遠超過電腦的負荷量，因此在 (14)中必須
對使用的項數做一些省略，稱為 limited-CI or truncated CI，最常用的 CI-SD 方
法就是只使用 (14) 中的前三項而忽略其他 highly excited configurations，更精確
但極為費時的CI-SDTQ 則是使用 (14) 中的前五項。單純的CI-SD 方法現在很
少使用在化學反應系統上，因為通常truncated CI 所得到的相對能量並不可靠，
一種較新的 quadratic CI 方法，如QCISD 或 QCISD(T) 方法則可以大幅改善相
對能量的問題。
CI 方法也可以很簡單的延伸到激發態的計算，比如說使用 (14) 中的前二
6
項並要求波函數與 HF 波函數 orthogonal，則可求得電子激發態的近似解，此
方法稱為 CIS (CI Singles)。若將 (14) 式中的第三項貢獻以 perturbation theory
的方式加入則可得到所謂的 CIS(D) 理論。 激發態的計算會在後面的章節做進
一步的介紹。
5. Coupled-Cluster Method (CC)
Coupled cluster 是早在1960 及1970 年代間發展出來的準確量子化學方法，
由於它的計算量很大，一直到了最近十年來才被廣泛的使用在一般的化學系統的
計算。CC 方法的基本方程式為
HF
Ψ = eTˆΨ (15)
= + + + +L
3 !
ˆ2 !
ˆˆ1
ˆT 2 T 3 eT T (16)
TˆTˆTˆTˆTˆn = 1 + 2 + 3 +L+ (17)
其中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
ˆˆ(18)
方程式 (15) 可被證明是得到完整的薛丁格方程式解的另一種寫法，而各項係數
t 則須由滿足 (15) 式的前提下經由複雜的推演而導出的一系列龐大的非線性方
程組中求反覆求解。完整的 CC 理論對於一般化學分子而言計算量過於龐大，
因此必須做一些簡化。通常的簡化方法是只保留 (17)式中的頭幾項而忽略高次
激發項。比如說，若只用 T2 則稱為 CCD 理論，若保留了 T1 及T2，則稱為 CCSD
理論。其實T3 項對於達到高準確度而言非常重要，只是若同時保留T1，T2，及
T3 則計算量通常會變的難以負荷，因此一種折衷的辦法是不直接使用T3 而是利
用 perturbation theory 等方法去估計T3 的貢獻，如現在常被用來當作高準確度
量化計算之標準方法的 CCSD(T) 理論就是屬於這種類型。一般認為 QCISD(T)
7
及 CCSD(T) 方法是當今常用的量化計算理論中最可靠的二種方法，不過它們必
須搭配適當的基底函數 (如 aug-cc-pVTZ) 才能完全發揮功能。
6. Density Functional Theory (DFT)
以薛丁格方程式求波函數是一項艱鉅的工作，在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 (19)
hi i us (ri )
2
ˆKS = −1 ∇2 + (20)
8
其中 us 是上述所謂的 external potential，θKS 稱為 density orbital 或
Kohn-Sham (KS) orbital。由此，我們可以證明系統的電子密度可以表示成
Σ=
=
n
i
i
1
KS 2 ρ θ (21)
而此虛擬系統中電子的動能可以寫成
Σ=
= −∇
n
i
Ks i i
1
2 KS
1
KS(1) | | (1)
2
1 θ θ (22)
電子與原子核間的位能可寫成
1
1
(r1) dr
r
UNe = −ΣZ ∫
α α
α
ρ
(23)
其中 α 代表不同的原子核，Zα 代表核電荷，r1α 代表空間中任一點到原子核 α
的距離，而電子之間的古典靜電位能可以寫成
1 2
12
(r1) (r2 ) r r
2
1 d d
r
Uee = ∫ ∫
ρ ρ
(24)
其中的1/2 是避免重複計算電荷間的作用力。 以上這些能量的計算都很直接，
只要有了上述的 Kohn-Sham orbitals 以及密度分佈函數 ρ(r) 就可以很容易的
得到。然而，式 (22) 與 (24) 有一些其他的問題。(22)是虛擬系統中的動能，真
實系統的動能必然有些不同；(23) 是古典靜電能量，並不包含量子力學中的電
子間的 exchange 及 correlation energy。我們若將這些校正量合稱為
exchange-correlation energy functional, Exc，則系統的總能量可以寫成
E =UNe + Ks +Uee + Exc (25)
Hohenberg and Kohn 也證明了如果我們知道真正的density functional，則真正的
電子密度會最小化系統的能量，這稱為 Hohenberg-Kohn variational theorem。因
此，我們可以藉著改變電子密度分佈函數來最小化式 (25) 所得到的能量，以求
得系統的電子密度與能量。由式 (21) 我們知道電子密度可由 KS orbitals 求
出，因此我們可以藉由改變 KS orbitals 來改變電子密度分佈函數。最佳化 (25)
所得到的能量我們可以得到一個類似Hartree-Fock Equation 的方程式稱為
Kohn-Sham Equation，可用類似 HF 理論中的 SCF 方法求解。如同 molecular
orbital，KS orbital 可以寫成是atomic orbital 或基底函數的線性組合，因此求解
Kohn-Sham Equation 也等於求構成KS orbital 最佳的基底函數線性組合。到目前
9
為止 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 (26)
而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 則在化學計算上的準確度可以大幅提升，此種方
法稱為 Hybrid DFT，其中最有名的莫過於在1993-1994 年所發展出的 B3LYP 方
法
LSDA LYP LSDA
0
HF B88
0
B3LYP (1 ) (1 ) Exc = a Ex + axEx + −a −ax Ex + acEc + −ac Ec
(27)
由與實驗值對照所得之最佳化係數為 a0 = 0.20, ax = 0.72, ac = 0.81。B3LYP 方法
是近十餘年來最受歡迎的 DFT 方法，因為幾乎在各種型態分子中的不同性質都
可以由此方法得到可靠的結果，許多人甚至將DFT 與 B3LYP 劃上等號。 然
而，最近幾年來發展出了一些新的 hybrid functionals 包括 mPW1PW91, B1B95,
B98, BMK, M05, M06, M08 等其準確度在很多方面已經可以超過 B3LYP。
Hybrid DFT 的計算量大約與系統大小的四次方成正比，與 Hartree-Fock 方法類
似，但現代的hybrid DFT 的準確度可以與MP2 方法相當甚至更好，一般認為
DFT 在未來的發展上應該還有很大的改進的空間，最近的研究也顯示，DFT 可
以與wavefunction 的方法相結合而得到很準確又有效率的理論方法。
胡維平
國立中正大學
化學暨生物化學系
?Copyright 2012
No comments:
Post a Comment