第五章 結構最佳化計算

第五章 結構最佳化計算

在位能曲面上的能量最低點我們稱為 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

計算所得的結果為：

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 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) CH₄

甲烷是最簡單的碳氫化合物，由於地球上生物圈的氣體自然循環以及石化燃料的 開採與使用，甲烷的在大氣中有相當的含量，也是造成溫室效應的一種重要氣體。甲 烷雖然有五個原子，但因為其高對稱性 (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

以 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) 光譜的譜線位置以及強度，我們在之後會對頻率的計算有更 進一步的說明。

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 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”(C₂) 結構，計 算所得之能量為MP2= 277.5807643，這二種結構的振動也都沒有虛頻。因此 C2 結構 比C_2h 結構在 MP2/6-31+G** 理論下約穩定 0.4 kcal/mol。因此對於較大的分子我們常 需要輸入各種可能的低能量結構來找出重要的能量最低點。

(F) 檢視電子密度及分子軌域

我們可以將計算得到的電子密度及分子軌域以圖形介面畫出來，在上例中，我們 產生了一個 phenol.chk 的 checkpoint file，其中包含了分子軌域與電子密度的資訊，

我們可以使用 formchk 的程式將其轉變成 formatted checkpoint file phenol.fchk 再將其 以文字方式傳到安裝 GaussView 的電腦上就可以各種不同的方式來呈現密度及軌域的 圖形，如圖三及圖四。

圖三 Phenol 分子的等電子密度表面 (0.0004 au) 顏色代表在表面上的不同的靜電位能。

圖四 Phenol 分子的 HOMO (左) 與 LUMO (右) 圖形

這些圖形雖然通常並不代表明確的物理性質，但在許多定性上的討論中仍有其參 考價值，或者至少穿插在報告或文章中看起來很漂亮。圖三中，氧原子附近帶有明顯 的負電位，而OH 基上的氫原子帶有較強的正電位，這與基本的化學觀念吻合。圖四 中，HOMO (highest occupied molecular orbital) 是一個  bonding orbital，而 LUMO (lowest unoccupied molecular orbital) 是一個  antibonding orbital，與預期的結果吻合；

注意在這二個軌域中，分子平面都是一個結面，而且軌域的數值變號對稱，這也是一 般  軌域的定義。

(G) 過渡態計算---Isomerization of Hydrogen cyanide

理論上我們也可以利用量子化學計算來估計化學反應的速率常數，因為課程時間 的關係，我們僅以一個很簡單的系統做能量障礙計算的說明：

CNH  HCN (R1)

我們使用 MP2/6-31+G(d,p) 計算方法得到：

E(HCN) = 93.17212 h E(CNH) = 93.14306 h

因此，Erxn = 0.02906 h = 18.2 kcal/mol，HCN 遠較 CNH 穩定。接下來我們嘗試計 算過渡態 (Transition State)：

%mem=200MW

%NProcShared=4

#MP2/6-31+G** OPT(CalcFC,TS,Noeigentest) FREQ HCN-->CNH TS

0 1 H C,1,RCH

N,2,RCN,1,AHCN Variables:

RCH=1.17 RCN=1.20 AHCN=74.0

在route section 中 OPT(CalcFC,TS,Noeigentest) 代表在尋找過渡態 (TS) 之前先求力常 數 (force constants)，這通常是在計算 TS 結構時必要的步驟。Noeigentest 告訴程式不 要因為估計出之eigenvalue 的符號或數目不對而過早放棄。TS 的起始結構的決定通常 必須由對反應中原子的移動的了解去估計，在這裏我們使用一個三角形的結構，並增

加 HCN 中的 CH 鍵長。計算結果是：

Final structure in terms of initial Z-matrix:

H

C,1,RCH

N,2,RCN,1,AHCN Variables:

RCH=1.16735021 RCN=1.1965236 AHCN=74.08555607

1 2 3 A' A' A'

Frequencies -- -1294.9335 2102.3079 2780.2898 Red. masses -- 1.1698 11.4434 1.0619 Frc consts -- 1.1557 29.7987 4.8363 IR Inten -- 205.5094 3.4266 202.2854 Atom AN X Y Z X Y Z X Y Z

1 1 -0.14 0.98 0.00 0.02 0.35 0.00 0.99 0.09 0.00 2 6 -0.07 -0.06 0.00 -0.03 0.70 0.00 -0.07 -0.01 0.00 3 7 0.07 -0.02 0.00 0.02 -0.62 0.00 -0.01 0.00 0.00

E(TS) = 93.08701 h

我們看到果然只有一個虛頻 (mode 1)，而且是代表著氫原子的移動。計算所得之能量 障礙 (Born-Oppenheimer energy) 為：

V^ = –93.08701 h – (–93.14306 h) = 0.05605 h = 35.17 kcal/mol

利用能量障礙以及TS 和反應物的結構及振動頻率等性質可以粗略估計出化學反應的 速率常數。

(H) 反應路徑計算---Isomerization of Hydrogen cyanide

尋找到過渡態後常需要進一步確認所找到的 TS 真的是所要研究的反應的過渡態，常

用的方法是利用 Gaussian 中 Intrinsic Reaction Coordinate (IRC) 功能的計算，沿著反應 路徑一步步搜尋在 IRC 上的分子結構及能量，看結構是不是逐漸收斂到反應物與產物，

能量變化是否與預期的反應吻合。 以前面的 HNC 反應為例：

%chk=HCN_IRC

%mem=120MW

HNC --> HCN IRC 0 1

H C 1 RCH

N 2 RCN 1 AHCN Variables:

RCH 1.16735 RCN 1.19652 AHCN 74.08556

或者直接利用前次的TS 計算結果：

%chk=HCN_TS

%mem=120MW

#MP2/6-31+G(d,p) Geom=Check, Guess=Check, IRC(RCFC,StepSize=5,MaxPoints= 15)

HNC --> HCN IRC 0 1

注意一定要使用 CalcFC 或 RCFC　取得 TS 的力常數才能成功的計算，也必須輸入正 確的 TS 結構 (CalcFC) 或者由入前次 TS 計算的 chk 檔讀入 (RCFC)。StepSize=5 代表 在 IRC 上每 0.05 bohr 計算一次結構，MaxPoints= 15 代表在 TS 前後各計算 15 點。

計算結束後可由 GaussView 軟體觀看 IRC 上結構及能量的變化，如上圖，一般定義 TS 的反應座標為零。

胡維平

國立中正大學 化學暨生物化學系

© Copyright 2018