MIT max, 为什么在我们意识体验所包含的信息内容似乎远大于37个比特? 化學鍵結: 2. 在最穩定鍵長的附近其實這二種能量貢獻都遠大於鍵能，所以化學鍵結其實是在量子化學中能量微妙平衡的結果。所以化學鍵結作为有效质量（老爱质能关系，有效质量=结合能除以光速平方），你化學鍵結有效质量太脆弱了。
消的結果，如下圖，在最穩定鍵長的附近其實這二種能量貢獻都遠大於鍵能，所. 以化學鍵結其實是在量子化學中能量微妙平衡的結果。由於電子的能量與me/h2.
分子性質與結構最佳化
由於原子可以互相鍵結形成各式各樣的化學分子，宇宙中才能有各式各樣豐 富的化學反應，生命體的存在更是仰賴複雜的化學環境與有機分子。但原子為何 會形成分子？這是以傳統化學觀念很不容易回答的問題，特別像是 H 2 分子， 為什麼二個相同的原子間可以有這麼強的化學鍵結？自從量子力學發展以來，化 學家便開始嘗試藉著解薛丁格方程式去了解化學分子鍵結的原理。
分子系統的薛丁格方程式通常非常複雜，因為其 Hamiltonian 包含了許多 項：
NeeeeNNN VVKVHK ˆˆˆˆˆˆ ++++= (1) KN 為原子核之動能，Ke 為電子之動能，VNN 為原子核間的靜電位能，VNe 為 原子核與電子間的靜電位能，Vee 為電子間的靜電位能。若要同時考慮原子核及 電子的運動將會使得問題非常複雜，一般在處理分子系統時我們都使用所謂的 Born-Oppenheimer approximation，這個假設的重點是：由於原子核比電子的質量 大很多，以電子的觀點而言，原子核幾乎是靜止的，而原子核的任何運動電子可 以在瞬間同時配合調整，因此我們可以將電子的運動與原子核的運動分開處理， 也就是說我們可以在固定的原子核位置下求電子的薛丁格方程式之解：
eleeleeleele Neeeelee ψEψH VVHK ˆ ˆˆˆˆ = ++=
(2)
因為此時 KN 可先設為零而 VNN為一常數。 解出 (2) 式後， 將 Eele 加上 VNN
U (R) = Eele (R) + VNN (R) (3)
則為在固定原子核位置下的總能量，或稱為 Born-Oppenheimer energy，其中 R 代 表原子核的座標。當然，原子核也會運動，如分子的振動、轉動等，分子轉動可 看成三度空間中的自由旋轉， 而分子的振動可看成是在由上述總能量所形成之 有效位能場下的運動。 U (R)常被稱為位能函數 potential energy function 或 potential energy surface (PES)。一般所談的分子結構可定義為 PES 上的能量最低 點所對應的原子核相對位置。
2
(1) H2+ (hydrogen molecular ion) 最簡單的化學分子是 H2+，此系統只有一個電子，在 Born-Oppenheimer 假 設下，此系統的波函數是可以用很複雜的數學方法解出來的，而得到的鍵長為 1.06 Å，鍵能為 2.8 eV。然而，這些數學方法並不適用於其他更複雜的系統，所 以我們只介紹由研究 H2+所得到的重要結果而不去深究其薛丁格方程式的解法。 由量子力學我們可以推論出線型分子的電子波函數同時也是分子軸方向角動量 的 eigenfunction：
hL 2, 1, 0, , ˆ ±±== mψmψL z (4)
由 H2+ 的解我們知道能量與 m2 有關， 因此若不考慮電子自旋，當 m ≠ 0，H2+ 的波函數是 doubly degenerate。仿照氫原子，我們將每個單電子波函數稱為一個 molecular orbital，而 m 值為 0, ±1, ±2 的軌域分別稱為 σ, π, δ 分子軌域。 在多 電子系統中，線性分子也可依照所帶的軸向角動量決定 term symbols (Σ, Π, Δ 等)。
由於此系統只有一個電子，所以 HF 方法已經是最佳理論。 下圖為使用 aug-cc-pVTZ basis set 計算所得的位能曲線
H2+ PES (HF/aug-cc-pVTZ)
-80
-60
-40
-20
0
20
40
60
0.5 1 1.5 2 2.5 3
R(H-H) (A)
E − E(H) (kcal/mol)
圖一 H2+ 的位能曲線
3
Gaussian 計算的 input file 如下:
%mem=200MW #HF/aug-cc-pVTZ SCAN
H2+ scan
1 2 H H 1 R
R=0.5 25 0.1
圖一中的縱軸為 H2+ 基態 (Σg) 總能量减掉一個基態氫原子的能量，由圖中 所預測的最穩定鍵長約為 1.06 Å，鍵能約為 64 kcal/mol (2.8 eV)，與實驗幾乎完 全吻合。圖一看似簡單，但卻隱藏著化學鍵結的一個重要概念，那就是原子間的 位能是原子核間的排斥力 (VNN) 與電子的動能與位能 (Eele = Ek + ENe) 互相抵 消的結果，如下圖，在最穩定鍵長的附近其實這二種能量貢獻都遠大於鍵能，所 以化學鍵結其實是在量子化學中能量微妙平衡的結果。由於電子的能量與 me/h2 成正比
2
0
2
4
4
1
ε
h me hartree e = (5)
H2+ PES (HF/aug-cc-pVTZ)
-1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2
0.5 1 1.5 2 2.5 3
R(H-H) (A)
E (hartree)
E(total) - E(H) V(NN) E(ele) - E(H)
圖二 H2+ 的位能曲線中原子核間的排斥力與電子能量的貢獻
4
如果蒲郎克常數再大二倍或者電子的質量再小二倍，穩定的化學鍵結可能就不會 存在了
(2) H2 (hydrogen molecule)
H2 是最簡單的中性分子，實驗上測得之 Re = 0.741 Å, De = 4.75 eV。雖然這 個系統只有二個電子，但由於電子的排斥力項使得我們無法得到薛丁格方程式的 解析解；，然而以現在電腦的能力我們其實可以很容易計算出在實驗誤差範圍內 H2 的任何性質。下圖為使用 CISD/aug-cc-pVTZ 方法計算所得的位能曲線
H2 PES (CISD/aug-cc-pVTZ)
-120
-100
-80
-60
-40
-20
0
0.5 1 1.5 2 2.5 3
R(H-H) (A)
E − 2 * E(H) (kcal/mol)
圖三 H2 的位能曲線
Gaussian 計算的 input file 如下:
%NProcShared=4 %mem=200MW #CISD(MaxCyc=200)/aug-cc-pVTZ SCAN
H2 scan
0 1 H H 1 R
5
R=0.5 25 0.1
(MaxCyc=200 的目的是讓 CISD 的計算有足夠的 cycle 數得到收斂的結果) 由於系統中只有二個電子，CISD 就等於 FCI 理論。圖二中的縱軸為 H2 基 態 (Σg) 總能量减掉二個基態氫原子的能量，也就是氫分子的鍵能。由圖中所預 測的最穩定鍵長約為 0.75 Å，鍵能約為 108 kcal/mol (4.7 eV)，與實驗值非常接 近。從圖三中我們觀察到幾件事：第一，在二個電子的系統中以非常簡單的計算 就可以達到非常高的精確度；第二，增加一個電子明顯的增加了電子與原子核間 的作用力，使得鍵能明顯的提升，鍵長明顯的縮短；第三，FCI 是一個 size-consistent 的理論，在鍵長很大時，總能量會等於二個氫原子之和。
圖三如圖二將鍵能的貢獻分成原子核間的排斥力 (VNN) 與電子的動能與位 能 (Eele = Ek + ENe + Eee)，由於電子的非定域 (delocalization) 性質， 此處 Eee 的貢獻遠小於 VNN 及 ENe。 就如同前述一個電子的系統，在最穩定鍵長的附近 VNN 和 Eele這二種能量貢獻都遠大於鍵能，這種情況在重原子的系統中更加明 顯。。
H2 PES (CISD/aug-cc-pVTZ)
-1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2
0.5 1 1.5 2 2.5 3
R(H-H) (A)
E (hartree)
E(total) - 2 * E(H) V(NN) E(ele) - 2 * E(H)
圖四 H2 的位能曲線中原子核間的排斥力與電子能量的貢獻
6
(3) 結構最佳化 (Geometry Optimization)
在位能曲面上的能量最低點我們稱為 energy minimum 或者是分子的穩定結 構，由前面二個例子我們看到可以藉著改變鍵長來找到能量的最低點，但這種方 法在多原子系統中 (包含很多鍵長與鍵角) 非常繁瑣，不易使用。如果理論方法 除了系統的能量外還能夠很有效率的提供 energy gradients 也就是能量梯度 (能 量對分子座標的一次微分)，那就會使得在多維空間中尋找能量最低點的工作 (結構最佳化) 容易很多。目前大部分的量化計算軟體中針對 HF, MP2, DFT 等 理論方法都可以提供所謂的 analytical gradients，也就是在計算完能量後很快速 的一次求出所有能量 gradients 的分量，如此可以配合很有效率的演算法 (如 Berny Algorithm) 在很少的步數下找到能量最低點。因此，MP2 及 DFT 方法常 被用來預測分子的穩定結構。然而目前在大部分的軟體中仍無法提供高準確度方 法如 QCISD(T) 或 CCSD(T) 等之 analytical gradients，因此使用這些方法做多 原子分子的結構最佳化仍然非常困難。以下以幾個簡單的例子來說明結構最佳化 的計算：
(A) H3+
H3+ 是星際空間中非常重要的離子，最主要的產生方式是：
H2+ + H2 → H3+ + H (6)
H3+ 是正三角形的結構，理論上我們也可以利用 scan 的方法做出如圖一或圖三 的位能曲線，然後求出能量最低點；但我們在這裡直接使用結構最佳化的方法來 做會更有效率：
#CISD/aug-cc-pVTZ OPT
H3+ Geometry Optimization (D3h)
1 1 H H 1 R H 2 R 1 60.
R=1.0
計算所得的結果為：
7
Final structure in terms of initial Z-matrix: H H,1,R H,2,R,1,60. Variables: R=0.87491145
HF=-1.2996582 CISD=-1.3417413
鍵長 (0.875 Å) 介於H2 和 H2+ 之間，而總鍵能為 −1.341 + 1.000 = 0.341 h = 214.0 kcal/mol，與實驗值 ~216 kcal/mol 非常接近，計算之平均 H−H 鍵能約為 71 kcal/mol，也是介於 H2 和 H2+ 之間。配合 analytical gradients 使用結構最佳 化的方法通常可以很快的將鍵長收斂到 0.001 Å 以下。此計算中的 correlation energy 佔總能量的 3%。
(B) H2O
水在地球上無所不在，但許多水的性質至今仍難以捉模。 氣態單一水分子 的性質可以使用量化的方法準確的模擬，比如說下例中我們計算水的結構，能 量，及其 dipole moment：
#MP2/aug-cc-pVTZ OPT Density=MP2
H2O (C2v)
0 1 H O 1 R H 2 R 1 A
R=0.95 A=104.5
Density=MP2 指定使用 MP2 的電子密度做進一步分析，主要結果如下：
Final structure in terms of initial Z-matrix: H O,1,R H,2,R,1,A Variables: R=0.96134417
8
A=104.13279222 HF=-76.0602858 MP2=-76.3289923
Mulliken atomic charges: 1 H 0.214777 2 O -0.429555 3 H 0.214777
Dipole moment (field-independent basis, Debye): X= 0.0000 Y= 0.0000 Z= -1.8591 Tot= 1.8591
實驗上得到的鍵長是 0.958 Å，鍵角是 104.5°，偶極距為 1.85 D。由於使用的理 論很簡單，能得到與實驗非常接近的結構以及準確的偶極距讓人感到相當的欣 慰。由於系統包含 10 個電子，已經沒有任何實用量化計算理論可以完整的描述 此系統的波函數；在此我們使用簡單的 MP2 理論配合還算大的基底函數 aug-cc-pVTZ 來計算分子的性質。注意此處的 correlation energy 只佔了總能量 的 0.3%。
(C) CH4
甲烷是最簡單的碳氫化合物，由於地球上生物圈的氣體自然循環以及石化燃 料的開採與使用，甲烷的在大氣中有相當的含量，也是造成溫室效應的一種重要 氣體。甲烷雖然有五個原子，但因為其高對稱性 (Td)，真正的結構參數只有 C−H 鍵長，但我們在輸入結構時要注意是否真的指定了我們想要輸入的 point group。比如說，使用 internal coordinates:
H C 1 R H 2 R 1 A H 2 R 1 A 3 120. H 2 R 1 A 3 -120.
R=1.09 A=109.4712
可以確定是 Td 結構，但如果設 A=109.5 則計算程式會將結構看成是 C3v。 請確認 output file 中的對稱性顯示
Stoichiometry CH4 Framework group TD[O(C),4C3(H)] Deg. of freedom 1 Full point group TD NOp 24 Largest Abelian subgroup D2 NOp 4 Largest concise Abelian subgroup D2 NOp 4
9
以 MP2/aug-cc-pVTZ 計算所得的最佳鍵長為 1.086 Å, 與實驗值 1.087 Å 幾乎 完全吻合。
(D) Phenol
對於較大的化學分子，使用 internal coordinates 建立分子座標會非常不方 便，通常我們會先使用包含圖形介面的化學軟體來畫出分子結構，再將其產生的 xyz 座標複製到 Gaussian 的 input file 中，以下我們以 phenol 分子 (C6H5OH) 為例，此分子包含了 13 個原子，用 internal coordinates 不是不可能，但太耗時 間也容易寫錯，我們於是先用 GaussView 軟體建好對稱性為 Cs 的結構，並將 其存成 xyz 格式的 gjf 檔，然後將結構選取貼至以下的 input file 中
%mem=400MW %chk=phenol %NProcShared=4 #B3LYP/6-31+G(d,p) OPT FREQ
Phenol Cs structure, and frequency calculation
0 1 C 0.06247519 -1.85905183 0.00000000 C 1.23928658 -1.10966130 0.00000000 C 1.17893066 0.28374414 0.00000000 C -0.05867147 0.92838157 0.00000000 C -1.23515866 0.17910532 0.00000000 C -1.17461093 -1.21471803 0.00000000 H 0.11027958 -2.95762252 0.00000000 H 2.21442367 -1.61794178 0.00000000 H 2.10646334 0.87448905 0.00000000 H -2.21068884 0.68686021 0.00000000 H -2.10222052 -1.80520041 0.00000000 O -0.12022972 2.35705594 0.00000000 H 0.77007233 2.71616879 0.00000000 其中 FREQ 代表除了結構最佳化外也計算分子振動頻率。計算頻率在此有二個 目的，第一是確認計算所得之結構確實為能量最低點，條件是所有的頻率皆為正 值；第二是預測此分子之紅外線 (IR) 光譜的譜線位置以及強度，我們在之後會 對頻率的計算有更進一步的說明。
10
1 2 3 A" A" A' Frequencies -- 229.1341 334.3462 405.5063 Red. masses -- 4.3526 1.1115 3.9446 Frc consts -- 0.1346 0.0732 0.3822 IR Inten -- 1.0012 114.8362 10.6479 計算結果顯示最小的頻率是 229 cm−1，並無虛頻，因此我們所輸入的 Cs 結構 的確是一個能量最低點。目前量化計算程式中大多對於 HF, MP2, DFT 等方法有 提供 analytical frequency (hessians)，可以很有效率的計算振動頻率與強度，但對 於高階的方法如 QCISD(T)，CCSD(T) 等，頻率的計算仍只能用數值方法求得， 對於較大的分子 ( > 4 個原子) 通常並不太實際。頻率的計算會同時得到振動模 式的 eigenvectors，可以藉由圖形介面來模擬每一種振動中原子移動的情形。另 外要注意的是頻率計算一定要使用與結構最佳化完全一樣的理論方法，否則得到 的結果通常沒有物理意義。
(E) 1,2-difluoro-ethane (CH2FCH2F)
在較大分子的結構最佳化計算中，我們要注意計算程式通常只會幫你找到距 離你輸入的起始結構最接近的能量最低點 (或稱穩定結構)，但此結構不見得是 此分子唯一的能量最低點，也很可能並不是能量最低的穩定結構。以下我們以 CH2FCH2F 分子為例：
%NProcShared=4 %mem=400MW #MP2/6-31+G** OPT FREQ
C2H4F2 (anti)
0 1 C C 1 R1 F 1 R2 2 A1 H 1 R3 2 A2 3 120. H 1 R4 2 A3 3 -120. F 2 R5 1 A4 3 180. H 2 R6 1 A5 3 60. H 2 R7 1 A6 3 -60.
R1=1.5 R2=1.3
11
R3=1.1 R4=1.1 R5=1.3 R6=1.1 R7=1.1 A1=110. A2=110. A3=110. A4=110. A5=110. A6=110.
此為一個“anti”(C2h) 結構，計算所得之能量為 MP2= −277.5800504。 若將第二個碳上的 F 與其中一個 H的二面角對調，則變為一個“gauch”(C2) 結構，計算所得之能量為 MP2= −277 .5807643，這二種結構的振動也都沒有虛 頻。因此 C2 結構比 C2h 結構在 MP2/6-31+G** 理論下約穩定 0.4 kcal/mol。因 此對於較大的分子我們常需要輸入各種可能的低能量結構來找出重要的能量最 低點。
(E) 檢視電子密度及分子軌域
我們可以將計算得到的電子密度及分子軌域以圖形介面畫出來，在上例中， 我們產生了一個 phenol.chk 的 checkpoint file，其中包含了分子軌域與電子密度 的資訊，我們可以使用 formchk 的程式將其轉變成 formatted checkpoint file phenol.fchk 再將其以文字方式傳到安裝 GaussView 的電腦上就可以各種不同 的方式來呈現密度及軌域的圖形，如圖四及圖五。
圖四 Phenol 分子的等電子密度表面 (0.0004 au) 顏色代表在表面上的不同的靜 電位能。
12
圖五 Phenol 分子的 HOMO (左) 與 LUMO (右) 圖形
這些圖形雖然通常並不代表明確的物理性質，但在許多定性上的討論中仍有 其參考價值，或者至少穿插在報告或文章中看起來很漂亮。圖四中，氧原子附近 帶有明顯的負電位，而 OH 基上的氫原子帶有較強的正電位，這與基本的化學 觀念吻合。圖五中，HOMO (highest occupied molecular orbital) 是一個 π bonding orbital，而 LUMO (lowest unoccupied molecular orbital) 是一個 π antibonding orbital，與預期的結果吻合；注意在這二個軌域中，分子平面都是一個結面，而 且軌域的數值變號對稱，這也是一般 π 軌域的定義。
胡維平 國立中正大學 化學暨生物化學系 © Copyright 2010
由於原子可以互相鍵結形成各式各樣的化學分子，宇宙中才能有各式各樣豐 富的化學反應，生命體的存在更是仰賴複雜的化學環境與有機分子。但原子為何 會形成分子？這是以傳統化學觀念很不容易回答的問題，特別像是 H 2 分子， 為什麼二個相同的原子間可以有這麼強的化學鍵結？自從量子力學發展以來，化 學家便開始嘗試藉著解薛丁格方程式去了解化學分子鍵結的原理。
分子系統的薛丁格方程式通常非常複雜，因為其 Hamiltonian 包含了許多 項：
NeeeeNNN VVKVHK ˆˆˆˆˆˆ ++++= (1) KN 為原子核之動能，Ke 為電子之動能，VNN 為原子核間的靜電位能，VNe 為 原子核與電子間的靜電位能，Vee 為電子間的靜電位能。若要同時考慮原子核及 電子的運動將會使得問題非常複雜，一般在處理分子系統時我們都使用所謂的 Born-Oppenheimer approximation，這個假設的重點是：由於原子核比電子的質量 大很多，以電子的觀點而言，原子核幾乎是靜止的，而原子核的任何運動電子可 以在瞬間同時配合調整，因此我們可以將電子的運動與原子核的運動分開處理， 也就是說我們可以在固定的原子核位置下求電子的薛丁格方程式之解：
eleeleeleele Neeeelee ψEψH VVHK ˆ ˆˆˆˆ = ++=
(2)
因為此時 KN 可先設為零而 VNN為一常數。 解出 (2) 式後， 將 Eele 加上 VNN
U (R) = Eele (R) + VNN (R) (3)
則為在固定原子核位置下的總能量，或稱為 Born-Oppenheimer energy，其中 R 代 表原子核的座標。當然，原子核也會運動，如分子的振動、轉動等，分子轉動可 看成三度空間中的自由旋轉， 而分子的振動可看成是在由上述總能量所形成之 有效位能場下的運動。 U (R)常被稱為位能函數 potential energy function 或 potential energy surface (PES)。一般所談的分子結構可定義為 PES 上的能量最低 點所對應的原子核相對位置。
2
(1) H2+ (hydrogen molecular ion) 最簡單的化學分子是 H2+，此系統只有一個電子，在 Born-Oppenheimer 假 設下，此系統的波函數是可以用很複雜的數學方法解出來的，而得到的鍵長為 1.06 Å，鍵能為 2.8 eV。然而，這些數學方法並不適用於其他更複雜的系統，所 以我們只介紹由研究 H2+所得到的重要結果而不去深究其薛丁格方程式的解法。 由量子力學我們可以推論出線型分子的電子波函數同時也是分子軸方向角動量 的 eigenfunction：
hL 2, 1, 0, , ˆ ±±== mψmψL z (4)
由 H2+ 的解我們知道能量與 m2 有關， 因此若不考慮電子自旋，當 m ≠ 0，H2+ 的波函數是 doubly degenerate。仿照氫原子，我們將每個單電子波函數稱為一個 molecular orbital，而 m 值為 0, ±1, ±2 的軌域分別稱為 σ, π, δ 分子軌域。 在多 電子系統中，線性分子也可依照所帶的軸向角動量決定 term symbols (Σ, Π, Δ 等)。
由於此系統只有一個電子，所以 HF 方法已經是最佳理論。 下圖為使用 aug-cc-pVTZ basis set 計算所得的位能曲線
H2+ PES (HF/aug-cc-pVTZ)
-80
-60
-40
-20
0
20
40
60
0.5 1 1.5 2 2.5 3
R(H-H) (A)
E − E(H) (kcal/mol)
圖一 H2+ 的位能曲線
3
Gaussian 計算的 input file 如下:
%mem=200MW #HF/aug-cc-pVTZ SCAN
H2+ scan
1 2 H H 1 R
R=0.5 25 0.1
圖一中的縱軸為 H2+ 基態 (Σg) 總能量减掉一個基態氫原子的能量，由圖中 所預測的最穩定鍵長約為 1.06 Å，鍵能約為 64 kcal/mol (2.8 eV)，與實驗幾乎完 全吻合。圖一看似簡單，但卻隱藏著化學鍵結的一個重要概念，那就是原子間的 位能是原子核間的排斥力 (VNN) 與電子的動能與位能 (Eele = Ek + ENe) 互相抵 消的結果，如下圖，在最穩定鍵長的附近其實這二種能量貢獻都遠大於鍵能，所 以化學鍵結其實是在量子化學中能量微妙平衡的結果。由於電子的能量與 me/h2 成正比
2
0
2
4
4
1
ε
h me hartree e = (5)
H2+ PES (HF/aug-cc-pVTZ)
-1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2
0.5 1 1.5 2 2.5 3
R(H-H) (A)
E (hartree)
E(total) - E(H) V(NN) E(ele) - E(H)
圖二 H2+ 的位能曲線中原子核間的排斥力與電子能量的貢獻
4
如果蒲郎克常數再大二倍或者電子的質量再小二倍，穩定的化學鍵結可能就不會 存在了
(2) H2 (hydrogen molecule)
H2 是最簡單的中性分子，實驗上測得之 Re = 0.741 Å, De = 4.75 eV。雖然這 個系統只有二個電子，但由於電子的排斥力項使得我們無法得到薛丁格方程式的 解析解；，然而以現在電腦的能力我們其實可以很容易計算出在實驗誤差範圍內 H2 的任何性質。下圖為使用 CISD/aug-cc-pVTZ 方法計算所得的位能曲線
H2 PES (CISD/aug-cc-pVTZ)
-120
-100
-80
-60
-40
-20
0
0.5 1 1.5 2 2.5 3
R(H-H) (A)
E − 2 * E(H) (kcal/mol)
圖三 H2 的位能曲線
Gaussian 計算的 input file 如下:
%NProcShared=4 %mem=200MW #CISD(MaxCyc=200)/aug-cc-pVTZ SCAN
H2 scan
0 1 H H 1 R
5
R=0.5 25 0.1
(MaxCyc=200 的目的是讓 CISD 的計算有足夠的 cycle 數得到收斂的結果) 由於系統中只有二個電子，CISD 就等於 FCI 理論。圖二中的縱軸為 H2 基 態 (Σg) 總能量减掉二個基態氫原子的能量，也就是氫分子的鍵能。由圖中所預 測的最穩定鍵長約為 0.75 Å，鍵能約為 108 kcal/mol (4.7 eV)，與實驗值非常接 近。從圖三中我們觀察到幾件事：第一，在二個電子的系統中以非常簡單的計算 就可以達到非常高的精確度；第二，增加一個電子明顯的增加了電子與原子核間 的作用力，使得鍵能明顯的提升，鍵長明顯的縮短；第三，FCI 是一個 size-consistent 的理論，在鍵長很大時，總能量會等於二個氫原子之和。
圖三如圖二將鍵能的貢獻分成原子核間的排斥力 (VNN) 與電子的動能與位 能 (Eele = Ek + ENe + Eee)，由於電子的非定域 (delocalization) 性質， 此處 Eee 的貢獻遠小於 VNN 及 ENe。 就如同前述一個電子的系統，在最穩定鍵長的附近 VNN 和 Eele這二種能量貢獻都遠大於鍵能，這種情況在重原子的系統中更加明 顯。。
H2 PES (CISD/aug-cc-pVTZ)
-1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 1.2
0.5 1 1.5 2 2.5 3
R(H-H) (A)
E (hartree)
E(total) - 2 * E(H) V(NN) E(ele) - 2 * E(H)
圖四 H2 的位能曲線中原子核間的排斥力與電子能量的貢獻
6
(3) 結構最佳化 (Geometry Optimization)
在位能曲面上的能量最低點我們稱為 energy minimum 或者是分子的穩定結 構，由前面二個例子我們看到可以藉著改變鍵長來找到能量的最低點，但這種方 法在多原子系統中 (包含很多鍵長與鍵角) 非常繁瑣，不易使用。如果理論方法 除了系統的能量外還能夠很有效率的提供 energy gradients 也就是能量梯度 (能 量對分子座標的一次微分)，那就會使得在多維空間中尋找能量最低點的工作 (結構最佳化) 容易很多。目前大部分的量化計算軟體中針對 HF, MP2, DFT 等 理論方法都可以提供所謂的 analytical gradients，也就是在計算完能量後很快速 的一次求出所有能量 gradients 的分量，如此可以配合很有效率的演算法 (如 Berny Algorithm) 在很少的步數下找到能量最低點。因此，MP2 及 DFT 方法常 被用來預測分子的穩定結構。然而目前在大部分的軟體中仍無法提供高準確度方 法如 QCISD(T) 或 CCSD(T) 等之 analytical gradients，因此使用這些方法做多 原子分子的結構最佳化仍然非常困難。以下以幾個簡單的例子來說明結構最佳化 的計算：
(A) H3+
H3+ 是星際空間中非常重要的離子，最主要的產生方式是：
H2+ + H2 → H3+ + H (6)
H3+ 是正三角形的結構，理論上我們也可以利用 scan 的方法做出如圖一或圖三 的位能曲線，然後求出能量最低點；但我們在這裡直接使用結構最佳化的方法來 做會更有效率：
#CISD/aug-cc-pVTZ OPT
H3+ Geometry Optimization (D3h)
1 1 H H 1 R H 2 R 1 60.
R=1.0
計算所得的結果為：
7
Final structure in terms of initial Z-matrix: H H,1,R H,2,R,1,60. Variables: R=0.87491145
HF=-1.2996582 CISD=-1.3417413
鍵長 (0.875 Å) 介於H2 和 H2+ 之間，而總鍵能為 −1.341 + 1.000 = 0.341 h = 214.0 kcal/mol，與實驗值 ~216 kcal/mol 非常接近，計算之平均 H−H 鍵能約為 71 kcal/mol，也是介於 H2 和 H2+ 之間。配合 analytical gradients 使用結構最佳 化的方法通常可以很快的將鍵長收斂到 0.001 Å 以下。此計算中的 correlation energy 佔總能量的 3%。
(B) H2O
水在地球上無所不在，但許多水的性質至今仍難以捉模。 氣態單一水分子 的性質可以使用量化的方法準確的模擬，比如說下例中我們計算水的結構，能 量，及其 dipole moment：
#MP2/aug-cc-pVTZ OPT Density=MP2
H2O (C2v)
0 1 H O 1 R H 2 R 1 A
R=0.95 A=104.5
Density=MP2 指定使用 MP2 的電子密度做進一步分析，主要結果如下：
Final structure in terms of initial Z-matrix: H O,1,R H,2,R,1,A Variables: R=0.96134417
8
A=104.13279222 HF=-76.0602858 MP2=-76.3289923
Mulliken atomic charges: 1 H 0.214777 2 O -0.429555 3 H 0.214777
Dipole moment (field-independent basis, Debye): X= 0.0000 Y= 0.0000 Z= -1.8591 Tot= 1.8591
實驗上得到的鍵長是 0.958 Å，鍵角是 104.5°，偶極距為 1.85 D。由於使用的理 論很簡單，能得到與實驗非常接近的結構以及準確的偶極距讓人感到相當的欣 慰。由於系統包含 10 個電子，已經沒有任何實用量化計算理論可以完整的描述 此系統的波函數；在此我們使用簡單的 MP2 理論配合還算大的基底函數 aug-cc-pVTZ 來計算分子的性質。注意此處的 correlation energy 只佔了總能量 的 0.3%。
(C) CH4
甲烷是最簡單的碳氫化合物，由於地球上生物圈的氣體自然循環以及石化燃 料的開採與使用，甲烷的在大氣中有相當的含量，也是造成溫室效應的一種重要 氣體。甲烷雖然有五個原子，但因為其高對稱性 (Td)，真正的結構參數只有 C−H 鍵長，但我們在輸入結構時要注意是否真的指定了我們想要輸入的 point group。比如說，使用 internal coordinates:
H C 1 R H 2 R 1 A H 2 R 1 A 3 120. H 2 R 1 A 3 -120.
R=1.09 A=109.4712
可以確定是 Td 結構，但如果設 A=109.5 則計算程式會將結構看成是 C3v。 請確認 output file 中的對稱性顯示
Stoichiometry CH4 Framework group TD[O(C),4C3(H)] Deg. of freedom 1 Full point group TD NOp 24 Largest Abelian subgroup D2 NOp 4 Largest concise Abelian subgroup D2 NOp 4
9
以 MP2/aug-cc-pVTZ 計算所得的最佳鍵長為 1.086 Å, 與實驗值 1.087 Å 幾乎 完全吻合。
(D) Phenol
對於較大的化學分子，使用 internal coordinates 建立分子座標會非常不方 便，通常我們會先使用包含圖形介面的化學軟體來畫出分子結構，再將其產生的 xyz 座標複製到 Gaussian 的 input file 中，以下我們以 phenol 分子 (C6H5OH) 為例，此分子包含了 13 個原子，用 internal coordinates 不是不可能，但太耗時 間也容易寫錯，我們於是先用 GaussView 軟體建好對稱性為 Cs 的結構，並將 其存成 xyz 格式的 gjf 檔，然後將結構選取貼至以下的 input file 中
%mem=400MW %chk=phenol %NProcShared=4 #B3LYP/6-31+G(d,p) OPT FREQ
Phenol Cs structure, and frequency calculation
0 1 C 0.06247519 -1.85905183 0.00000000 C 1.23928658 -1.10966130 0.00000000 C 1.17893066 0.28374414 0.00000000 C -0.05867147 0.92838157 0.00000000 C -1.23515866 0.17910532 0.00000000 C -1.17461093 -1.21471803 0.00000000 H 0.11027958 -2.95762252 0.00000000 H 2.21442367 -1.61794178 0.00000000 H 2.10646334 0.87448905 0.00000000 H -2.21068884 0.68686021 0.00000000 H -2.10222052 -1.80520041 0.00000000 O -0.12022972 2.35705594 0.00000000 H 0.77007233 2.71616879 0.00000000 其中 FREQ 代表除了結構最佳化外也計算分子振動頻率。計算頻率在此有二個 目的，第一是確認計算所得之結構確實為能量最低點，條件是所有的頻率皆為正 值；第二是預測此分子之紅外線 (IR) 光譜的譜線位置以及強度，我們在之後會 對頻率的計算有更進一步的說明。
10
1 2 3 A" A" A' Frequencies -- 229.1341 334.3462 405.5063 Red. masses -- 4.3526 1.1115 3.9446 Frc consts -- 0.1346 0.0732 0.3822 IR Inten -- 1.0012 114.8362 10.6479 計算結果顯示最小的頻率是 229 cm−1，並無虛頻，因此我們所輸入的 Cs 結構 的確是一個能量最低點。目前量化計算程式中大多對於 HF, MP2, DFT 等方法有 提供 analytical frequency (hessians)，可以很有效率的計算振動頻率與強度，但對 於高階的方法如 QCISD(T)，CCSD(T) 等，頻率的計算仍只能用數值方法求得， 對於較大的分子 ( > 4 個原子) 通常並不太實際。頻率的計算會同時得到振動模 式的 eigenvectors，可以藉由圖形介面來模擬每一種振動中原子移動的情形。另 外要注意的是頻率計算一定要使用與結構最佳化完全一樣的理論方法，否則得到 的結果通常沒有物理意義。
(E) 1,2-difluoro-ethane (CH2FCH2F)
在較大分子的結構最佳化計算中，我們要注意計算程式通常只會幫你找到距 離你輸入的起始結構最接近的能量最低點 (或稱穩定結構)，但此結構不見得是 此分子唯一的能量最低點，也很可能並不是能量最低的穩定結構。以下我們以 CH2FCH2F 分子為例：
%NProcShared=4 %mem=400MW #MP2/6-31+G** OPT FREQ
C2H4F2 (anti)
0 1 C C 1 R1 F 1 R2 2 A1 H 1 R3 2 A2 3 120. H 1 R4 2 A3 3 -120. F 2 R5 1 A4 3 180. H 2 R6 1 A5 3 60. H 2 R7 1 A6 3 -60.
R1=1.5 R2=1.3
11
R3=1.1 R4=1.1 R5=1.3 R6=1.1 R7=1.1 A1=110. A2=110. A3=110. A4=110. A5=110. A6=110.
此為一個“anti”(C2h) 結構，計算所得之能量為 MP2= −277.5800504。 若將第二個碳上的 F 與其中一個 H的二面角對調，則變為一個“gauch”(C2) 結構，計算所得之能量為 MP2= −277 .5807643，這二種結構的振動也都沒有虛 頻。因此 C2 結構比 C2h 結構在 MP2/6-31+G** 理論下約穩定 0.4 kcal/mol。因 此對於較大的分子我們常需要輸入各種可能的低能量結構來找出重要的能量最 低點。
(E) 檢視電子密度及分子軌域
我們可以將計算得到的電子密度及分子軌域以圖形介面畫出來，在上例中， 我們產生了一個 phenol.chk 的 checkpoint file，其中包含了分子軌域與電子密度 的資訊，我們可以使用 formchk 的程式將其轉變成 formatted checkpoint file phenol.fchk 再將其以文字方式傳到安裝 GaussView 的電腦上就可以各種不同 的方式來呈現密度及軌域的圖形，如圖四及圖五。
圖四 Phenol 分子的等電子密度表面 (0.0004 au) 顏色代表在表面上的不同的靜 電位能。
12
圖五 Phenol 分子的 HOMO (左) 與 LUMO (右) 圖形
這些圖形雖然通常並不代表明確的物理性質，但在許多定性上的討論中仍有 其參考價值，或者至少穿插在報告或文章中看起來很漂亮。圖四中，氧原子附近 帶有明顯的負電位，而 OH 基上的氫原子帶有較強的正電位，這與基本的化學 觀念吻合。圖五中，HOMO (highest occupied molecular orbital) 是一個 π bonding orbital，而 LUMO (lowest unoccupied molecular orbital) 是一個 π antibonding orbital，與預期的結果吻合；注意在這二個軌域中，分子平面都是一個結面，而 且軌域的數值變號對稱，這也是一般 π 軌域的定義。
胡維平 國立中正大學 化學暨生物化學系 © Copyright 2010
No comments:
Post a Comment