第七章 熱力學性質計算

第七章 熱力學性質計算

熱力學 (thermodynamics) 是探討巨觀世界中物質與能量之間的關係，也幾乎可以 說是用來描述自然界中現象最實用的一門科學。理論上，巨觀的現象可以用微觀世界中 分子的性質以及行為來解釋，但由於巨觀物體中包含非常大量的分子，因此實際上在巨 觀世界中觀察到的性質與現象並非由單獨分子的行為來決定，而是大量分子的行為之平 均效應。因此，我們可以藉由統計熱力學 (statistical thermodynamics) 的理論將分子的 性質與熱力學的性質聯繫起來；然而其中最困難的地方在於處理分子之間複雜的作用力，

所以凝態的性質並不容易以小分子的計算來準確模擬，但對於低壓氣態分子系統而言，

理想氣體 (ideal gas) 是很好的模型，由於在理想氣體中分子間沒有作用力，我們可以 很容易的利用統計方法從化學分子的量化計算結果推導出氣態分子的各種熱力學性質。

本章中關於熱力學性質計算的介紹以理想氣體為主。

首先我們要先來談一下能量的零點；一般而言能量都是相對的，在絕大部份的問題 中只有能量差是有意義的，因此在熱力學中並不需要有所謂的絕對能量。然而，在計算 化學中，訂出一個適用於大部分情況的能量零點在處理上及討論上是比較方便的，因此 一般是將所有的電子與原子核分離到無限遠並且是相對靜止時的能量狀態定為量化計 算中的能量零點。這種定義與前面所提到在 Born-Oppenheimer 假設下解薛丁格方程式 的能量零點是一致的，因此以下的討論中我們可以把分子中原子核運動的能量直接加到 Born-Oppenheimer energy 中，當作分子的總能量。

氣態分子系統的能量除了 Born-Oppenheimer energy 外主要來自於分子的各種運動，

如質心的直線運動 (translational motion)，分子的轉動 (rotation)，以及分子內原子間的 振動等 (vibration)。

(1) 質心的直線運動

分子在空間中質量中心的直線運動其實就是 particle in a 3-D box 的問題，能階為：











  







 2

2 2

2 2

2 2

8 c

n b n a n E m

E E

E_{trans x y z}  ^{x y z}

(1)

在巨觀系統中，分子可以活動的空間非常大，因此 a,b,c 都是很大的數，能階間的差異 非常小，幾乎可以看成是連續的能量分佈。若假設分子在各能階的數量是依照波茲曼分 佈 (Boltzmann distribution)，我們可以推導出分子的平均動能為

E̅_trans = ³₂kT = ³₂RT/mol (2)

此結果與古典力學所預測的一致。在 298.15 K 時，RT 能量為 0.592 kcal/mol，所以一 莫爾分子的 translational energy 為 0.889 kcal/mol。

(2) 分子的轉動

在 rigid rotor 的假設下，分子的旋轉能階為：

E_rot =_2I^Ja2

a+_2I^Jb2

b+_2I^Jc2

c (3)

I_a, I_b, I_c 為三個主軸 (principal axes) 上的轉動慣量，它們可以由最佳幾何結構輕易的求 出，J_a, J_b, J_c 則為對應的量子化角動量。由於分子轉動的能階差也非常小，除了在極低 溫下，轉動能量也可以看成是連續的，由波茲曼分佈我們也可以得到：

E̅_rot =³₂kT =³₂RT/mol (for nonlinear molecules) (4)

E̅_rot = kT = RT/mol (for linear molecules) (5)

所以在 298.15 K 時一莫爾非線性分子的 rotational energy 也是 0.889 kcal/mol。

線性分子軸上的轉動不會被激發所以通常不需考慮它的能量貢獻，因此轉動能量上只有 二個自由度需要考慮。

(3) 分子的振動

對於一個包含 n 個原子的非線性分子，原子核的運動包括了三個 translational motions，三個 rotational motions；但總共應有 3n 個自由度，所以剩下的 3n – 6 個 vibrations。分子的振動比較難處理，一般做法是利用 normal mode analysis 將這些振動 分解成 3n – 6 個獨立的 harmonic oscillators，而振動的總能量即為這些獨立振動的能量 和。分子的振動與直線運動和轉動有二點很大的差異；首先，分子的振動包含有振動零 點能 (vibrational zero-point energy, ZPE)，在簡諧振動假設下是

1/2 h，這是無法移除的一種能量；另外，振動的能階差 (h in harmonic approximation) 遠較另外二種運動來的大很多，在室溫下只有很少量的分子處於振動激發態，因此振動 能量不能用古典力學方法來處理，而需要用量子統計的方法。在 harmonic 假設下，振 動的能階為：

 2, 0,1, v

2 , v 1

,  



 



 

 _{i i i}

i

vib h

E  (6)

其中 i 指第 i 個振動模式。由統計熱力學，溫度為 T 時的平均振動能量是：

2 1 1

, /

 

 h kT

i vib

e h h

E 



而振動總能量則為：





^{3 6}

1 , n i

i vib

vib

E

E

在 0 K 下只有 ZPE 的貢獻，0 K 以上 (7) 式的第二項，也就是 vibrational thermal excitation energy，也會隨著溫度增高貢獻逐漸增大，但在室溫以下，絕大部分的分子都 是處於振動的基態，(7) 式右邊第二項的貢獻通常不大。 振動零點能可以儲存很大的 能量，比如說像甲烷這麼小的分子就包含了約 25 kcal/mol 的 ZPE，然而在室溫下，振 動激發能量只有約 0.2 kcal/mol。

(4) 內能與焓 (Internal Energy and Enthalpy)

內能 (U) 在熱力學中扮演著非常重要的角色，熱力學第一定律：

U = q + w (9)

就是內能，功與熱能的關係。 從分子層面來看，處於電子基態的理想氣體分子的內能 可寫成是：

vib rot

trans O

B E E E

E

U  _{. .}   (10)

等號右邊第一項是指由量化計算所得之 Born-Oppenheimer energy。

焓 (H) 也是在熱力學中經常用到的狀態函數，由定義：

H = U + PV

再配合理想氣體公式 PV = nRT，我們得到一莫爾的理想氣體分子的焓為：

RT E

E E

E PV U

H    B_.O_. trans rot  vib (12)

當然此時上式的每一項能量都要以每莫爾為單位，注意 Born-Oppenheimer energy 是不 隨溫度而改變的。當計算一個化學反應的反應熱時，我們可以個別計算產物及反應物的 焓，然後將產物焓的總和減去反應物焓的總和就可以得到 enthalpy of reaction (H_rxn)；

原則上H_rxn 會是溫度的函數，但通常其值對溫度的變化不會太敏感；從以上的討論也 不難看出，在 0 K 時，H_rxn 就等於U_rxn。

(5) 亂度 (Entropy)

在熱力學中，亂度是指能量在各種可能的能階中的分佈自由度，由統計熱力學我們 可以推出：

Q T k

S E

Q T k

ZPE E

T S U

Q T k

K U T S U

thermal O B

ln ) ln

(

) ln 0 ( ) (

. .





 

 

 



(13)

其中 Q 為系統的 canonical partition function，請參閱任何一本物化教科書，E_thermal 代 表溫度由 0 K 上升至 T 所增加的內能。在大部分情況下 Q 可以寫成是分子各種運動 的 partition function 之乘積：

vib rot trans

Q Q Q

Q

其中各種 partition function 都可以由各種運動的能階求得。由於 U 和 Q 都可以分解成 來自各種運動的貢獻，S 也可以寫成

vib rot trans

thermal S S S

T

S  E    (15)

一般而言，對相同的分子溫度增高會讓上式所有項都增大，體積加大也會讓 S_trans 增 加，較低的振動頻率也會對 S_vib 有比較大的貢獻。以一莫爾最簡單的單原子理想氣體 為例:

2 2 3 3 2

ln 2 2 ln

ln 3 2

5 2 ln 3











 











mk h R N

T R V R R S

Q T k

RT S

A trans



(16)

此公式稱為 Sackur-Tetrode equation，可以對鈍氣原子的亂度作非常準確的預測，比如 說由 (16) 式所得之 Ne 在 298 K 標準狀態下 (1 atm) 的亂度為 146.2 J mol^1 K^1，而 實驗值為 146.5 J mol^1 K^1。值得注意的是 (16) 式中第二項與體積有關，這代表莫爾 亂度為壓力的函數，因此就算是理想氣體，壓力不同，亂度也不同，壓力低，體積大，

亂度就高。

(6) Gibbs 自由能 (Gibbs Free Energy) 與平衡常數

Gibbs 自由能是一般討論化學反應的自發性及化學平衡最重要的函數，由定義：

G = H  TS (17)

對理想氣體而言 G 可由前述之 H 及 S 函數得到。熱力學第二定律告訴我們自發反應 的 G < 0，而平衡常數與G 有以下的關係：

Keq RT G  ln

 ^ (18)

在應用上，很多人對此公式是常有一些誤解，等號左邊是指所有反應物及產物都在標準 狀態下時的自由能變化，等號右邊的平衡常數是指在平衡時以相對於標準狀態的濃度比 例算出之無單位的濃度比值。若反應前後的分子數目不等，則要特別注意所得出之平衡 常數的單位，此單位會是使用標準濃度下的平衡常數單位，例如 2A  B 的反應，算 出之單位通常為 atm^1。 通常若反應物及產物的分子個數相等，H° 和 G° 差別不 大，但個數若不同，由於亂度的效應，在溫度較高時，分子個數較多的一側會有明顯較 低的 G 值而造成 H° 和 G° 很大的差異。理論上，平衡常數可由 (18) 式預測出，

但在實際運用上卻有些困難。首先，平衡常數對 G° 非常敏感，就算計算值只有 2 kcal/mol 的誤差 (這已經是目前很好的計算準確度)，室溫下 K_eq 值就會有 20 倍的誤差。

此外，在溶液中由於 solvation 的效應，使得準確的能量以及亂度的計算非常困難。有 關溶劑效應對於平衡的影響會在後面的章節中做進一步的討論。

(7) 莫爾生成熱、焓、與鍵能的關係

分子的莫爾生成熱在熱力學中是非常重要的資料，因為由黑斯定律我們知道任何化 學反應的反應熱可以由生成物及反應物的莫爾生成熱求得，而從熱力學的表格或資料庫 中可以查到許多常見的化學物質的莫爾生成熱。莫爾生成熱的定義是在標準狀態下由最 穩定的元素態反應成分子的過程中焓的變化。從計算化學的觀點來看，從能量零點出發 到形成氣態分子這包含了三個步驟：

A. 能量零點  組成分子的原子 H_atoms

=H

如果元素態是氣態，我們當然可以直接由第三步求得化合物分子與元素態的 H 値然後 相減得到H_f，然而我們也可以將 –[(B)+(C)] 的能量變化看成是化合物分子的總鍵能 (total atomization energy, TAE)，TAE(molecule) ，並將 (B) 的能量變化看成是元素態 分子的總鍵能 TAE(elements)，負號的來源是因為一般將鍵能定為正值。總鍵能的求法 非常容易，只須將計算所得之分子能量減掉原子能量 (步驟 A) 即可。這樣分解的好處 是：雖然 (C) 步驟中分子的種類幾乎可以說是無窮多，(A) 與 (B) 步驟的原子及元素 態的種類卻是非常有限。對於常見的 H, C, N, O , F, S, P, Cl 等元素 (B) 步驟已有非常 準確的實驗值，(A) 步驟也容易得到準確的計算值，所以當元素態是氣態時，分子的 H_f 可以由 TAE(molecule)之計算值與 TAE(elements)的實驗值而得到，此處之 TAE 代表 enthalpy 的變化。

H_f = TAE_calc (molecule)  TAE_expt(elements) (19)

當然，理論上 TAE(elements) 也可以計算出來，但因為在標準狀態下有些元素態是固體，

如碳，磷，及金屬等， (C) 步驟必須包含元素態的氣化能量，此能量不易計算，通常 只能由實驗得到。比如說對水分子而言，

H_f (H₂O) = TAE_calc (H₂O)  TAE_expt(H₂) + 1/2 TAE_expt(O₂)] (20)

但對甲烷分子而言：

H_f (CH₄) = TAE_calc (CH₄)  2 TAE_expt(H₂) + H_vap(graphite)] (21)

若上述之計算或理論值都是在一大氣壓下，則所得之結果為標準狀態下之值H°_f 。在 計算化學中 TAE 的預測也常被用來當作理論方法的測試，有時候我們也可以將準確估 計的振動零點能和熱能的貢獻從實驗測量的 TAE 中移除得到純粹只包含

Born-Oppenheimer energy 的 TAE_B.O.。使用 TAE_B.O. 來測試理論方法的好處是我們可 以單純針對電子薛丁格方程式的解做比較，而不受處理原子核運動 (分子振動及轉動) 所使用的假設而造成的不準度的影響。

(8) 計算實例的說明

以下我們以幾個例子來說明前述的熱力學性質的計算，首先我們以水為例：

Gaussian input file 如下：

%NProcShared=4

%chk=H2O

%mem=400MW

#MP2/aug-cc-pVTZ OPT Density=MP2 FREQ H2O (C2v)

0 1 H O,1,R H,2,R,1,A Variables:

R=0.96134417 A=104.13279222

我們使用 MP2 的理論配合 aug-cc-pVTZ 基底函數，這對水分子而言應該是一個很可 靠的計算。要得到熱力學資料，我們必須先做結構最佳化，然後再進行 frequency 的計 算。

計算完成後，首先我們應該確認這的確是一個結構最佳化收斂的結果：

Item Value Threshold Converged?

Maximum Force 0.000063 0.000450 YES RMS Force 0.000038 0.000300 YES Maximum Displacement 0.000219 0.001800 YES RMS Displacement 0.000220 0.001200 YES Predicted change in Energy=-1.361821D-08

Optimization completed.

-- Stationary point found.

看來結構已經收斂沒有問題，這的確是在 MP2/aug-cc-pVTZ 理論下的靜態點所做的頻 率計算。接下來我們看頻率計算的結果：

1 2 3 A1 A1 B2 Frequencies -- 1628.1074 3822.1217 3948.1174 Red. masses -- 1.0828 1.0450 1.0815 Frc consts -- 1.6911 8.9947 9.9324 IR Inten -- 71.7189 5.5610 75.5556 Atom AN X Y Z X Y Z X Y Z 1 1 0.00 0.43 0.56 0.00 0.59 -0.40 0.00 -0.56 0.43 2 8 0.00 0.00 -0.07 0.00 0.00 0.05 0.00 0.07 0.00 3 1 0.00 -0.43 0.56 0.00 -0.59 -0.40 0.00 -0.56 -0.43

我們得到三個振動頻率 (cm^1)，這是我們所預期的，因為 3N  6 = 3。第一個頻率是指 HOH bending，第二、三個頻率分別是 OH symmetric 和 antisymmetric stretching。

使用 MP2 理論配合簡諧振動假設通常會略為高估 (~10%) 振動頻率。

接下來是一些與熱化學相關的分子轉動與振動的資訊：

--- - Thermochemistry - ---

Temperature 298.150 Kelvin. Pressure 1.00000 Atm.

Atom 1 has atomic number 1 and mass 1.00783 Atom 2 has atomic number 8 and mass 15.99491 Atom 3 has atomic number 1 and mass 1.00783 Molecular mass: 18.01056 amu.

Principal axes and moments of inertia in atomic units:

1 2 3 Eigenvalues -- 2.23264 4.13829 6.37093 X 0.00000 0.00000 1.00000 Y 1.00000 0.00000 0.00000 Z 0.00000 1.00000 0.00000 This molecule is an asymmetric top.

Rotational symmetry number 2.

Rotational temperatures (Kelvin) 38.79429 20.92988 13.59516 Rotational constants (GHZ): 808.34280 436.10846 283.27757

Zero-point vibrational energy 56214.6 (Joules/Mol) 13.43561 (Kcal/Mol) Vibrational temperatures: 2342.48 5499.17 5680.45 (Kelvin)

首先，在 Gaussian 中熱力學計算內定的狀態是 298.15 K, 1 atm，這可由

Temperature 和 Pressure keywords 來改變。其次是每個原子的質量，注意此處使用的是 純的同位素質量而不是常用的平均原子量，若要變更成別種同位素可使用 Freq(ReadIso) keyword。比如說，若我們完成了 H₂O 的計算，我們可以用以下的 input file 來計算 D₂O 的熱力學性質。

%chk=H2O

#MP2/aug-cc-pVTZ Geom=Check FREQ(ReadFC,ReadIso) D2O (C2v)

0 1

298.15 1.0 2

16 2

接下來 output file 會列出分子的轉動慣量 (moments of inertia)，Gaussian 印出的是 沿著三個互相垂直的主軸 (principal axes) 的值 (I_A, I_B , I_C )，單位是 (amu) (bohr)2。 轉 動溫度與轉動常數的定義是：

k I h

r 8 ²

2

 

 (22)

8 ²

I B h





在簡諧運動的假設下計算出之 ZPE 為 13.4 kcal/mol，而振動溫度的定義是：

k h

v





 (24)

這些量對計算熱力學性質或分析及光譜性質有時候會比較方便。接著是熱力學函數的計 算結果：

Zero-point correction= 0.021411 (Hartree/Particle) Thermal correction to Energy= 0.024246

Thermal correction to Enthalpy= 0.025191 Thermal correction to Gibbs Free Energy= 0.003767

Sum of electronic and zero-point Energies= -76.307581 Sum of electronic and thermal Energies= -76.304746 Sum of electronic and thermal Enthalpies= -76.303802 Sum of electronic and thermal Free Energies= -76.325225

E (Thermal) CV S KCal/Mol Cal/Mol-Kelvin Cal/Mol-Kelvin Total 15.215 6.009 45.090 Electronic 0.000 0.000 0.000 Translational 0.889 2.981 34.608 Rotational 0.889 2.981 10.475 Vibrational 13.437 0.048 0.007

首先是將 ZPE 以 hartrees 來表示，第二行是 ZPE 加上 (13) 式中的 thermal energy，第 三行是第二行再加上 RT，第四行式在加上溫度乘上亂度 (T × S) 。第五到第八行則是 將計算出之 Born-Oppenheimer energy (MP2 = 76.3289923) 分別加到第一到四行的能 量而得到 U(0 K), U(298.15 K), H(298.15 K), 以及 G(298.15 K)。 接下來程式印出各種 運動對內能、定容比熱，以及亂度的貢獻。Translation 和 Rotation 對內能的貢獻就是 由 (2) 和 (4) 式所得到的 1.5 RT`; Vibration 對內能的貢獻就是由 (7) 和 (8) 式所得到 振動能量，注意此處的能量包含 ZPE，並不只是熱能。定容比熱 (Cv) 的定義是：

v

v 



 







  T

C U (25)

C

gas) 推出解析公式。實驗上 C_v (298.15 K) = 6.04 cal/mol-K 與此處的計算值 6.01 非常 接近! 亂度也可以由 (15) 式配合統計熱力學得到，由 output 可以看出，亂度最主要的 貢獻是來自於 translational motion 其次是 rotational motion，在室溫下 vibration 幾乎沒 有貢獻；實驗上 S (298.15 K) = 45.1 cal/mol-K 與此處的計算值幾乎完全吻合! 說實話，

此處 C_v和 S 計算的準確度真的有點讓人難以置信。 其實，C_v的準確度是可以預期的，

因為對理想氣體而言，translation 和 rotation 的內能和是 3RT，由 (25) 式得知它們對 比熱的貢獻是 3R = 5.96 cal/mol-K，vibration 的能階在室溫下被激發的量很少，所以對 比熱只有些微的貢獻，至於 S 的計算結果可以看成是量子力學與統計熱力學的成功的 驗證，沒有量子力學就不可能得到準確的絕對亂度值。

當一個化學反應的產物分子數目多於反應物時 (例如 dissociation, fragmentation reactions)，部分的 vibration modes 轉變成 translational 以及 rotation motions，因此亂 度會大幅增加，G_rxn的值在較高溫度下會遠低於 U_rxn 和 H_rxn。反過來說，若一個 化學反應的產物分子數目少於反應物時 (例如 association, dimerization, polymerization reactions)，G_rxn的值在較高溫度下會遠高於 U_rxn 和 H_rxn。

接下來我們來計算水分子的莫爾生成熱，由 (20)：

H°_f (H₂O) = TAE_calc (H₂O)  TAE_expt(H₂) + 1/2 TAE_expt(O₂)]

= TAE_B.O. (H₂O) + 2 E_thermal (H) + E_thermal (O)  ZPE(H₂O)  E_thermal (H₂O) + 2 RT –

[ 2 H°_f (H) + H°_f (O) ] (26)

2 RT 項為由U 轉換成 H 的校正項，因為由水分子變成三個原子，粒子數多了二。

TAE_calc (H₂O) 可由計算所得之 MP2/aug-cc-pVTZ 能量求得：

E(H) = 0.49982 h; E(O) = 74.95929 h; E(H₂O)= 76.32899 h TAE_B.O. (H₂O) = 0.42336 h = 232.2 kcal/mol

一個原子的 thermal energy 為 1.5 RT，由 output 中我們得到水分子的 ZPE + E_thermal =

15.2 kcal/mol，因此 (26) 可以化簡為：

H°_f (H₂O) = 232.2 kcal/mol + 4.5 RT  15.2 kcal/mol + 2 RT –

[ 2 H°_f (H) + H°_f (O) ] (27) 由 298.15 K 的實驗值：H°_f (H) = 52.1 kcal/mol; H°_f (O) = 59.6 kcal/mol

我們可以從 (27) 得到在 298.15 K 時 H°_f (H₂O) = 57.0 kcal/mol，而實驗值是 57.8 kcal/mol，誤差可以說是令人驚訝的小，因為通常在 TAE 計算時，MP2 會低估鍵能，

此處由於誤差的抵銷使得結果非常接近實驗值，一般而言，以單一 MP2 理論計算生成 熱會有數仟卡的誤差。通常我們會需要更高階的計算，例如 CCSD(T)/aug-cc-pVQZ 或 高精確度理論，如 G4，W1, MCG3-DFT，MLSE-DFT 等才能將鍵能算得很準確。

胡維平

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

© Copyright 2018