摘要：本教程介绍了焓、焓变、标准摩尔生成焓和标准摩尔反应焓的相关概念，并以顺式和反式1,3-丁二烯为例，详细介绍了如何利用Gaussian程序计算反应的标准摩尔反应焓。此外，以HCl分子为例，演示了如何计算分子的标准摩尔生成焓。最后，还演示了如何用原子化法计算生成焓。

作者：陈宇
日期：2020-03-13

基本概念

焓（enthalpy）

在热力学中，焓是用来表征系统的状态参数，当一个系统的状态不变时，系统有唯一的焓值，其通常用字母H表示，其定义为：

H = U + pV

其中U为体系总的内能，p为体系的压力，V为体系的体积。

焓变（enthalpy change）

上面我们已经阐述了焓的概念，当系统从一个状态变成另一个状态时，焓也会相应的从H(initial)变为H(final)，其焓值的改变量称为焓变，用字母ΔH表示：

ΔH = H(final) – H(initial)

标准摩尔生成焓（standard enthalpy of formation）

由于一个体系的绝对焓值很难测定，因此在实际应用中，我们只使用其相对焓值，对于一个分子而言，我们通常用标准摩尔生成焓来表征其焓值，用符号ΔfHmΘ表示，其定义为在标准状态下，由最稳定的单质生成该物质的焓变。其中标准态是指1 atm和298.15 K。

标准摩尔反应焓（standard molar enthalpy of reaction）

对于一个化学反应而言，如果它是在标准状态下进行，那么其反应焓变就称为标准摩尔反应焓变，用符号ΔrHmΘ表示。实际上标准摩尔生成焓就是一种特殊的标准摩尔反应焓，只不过其反应物必须是标准态下最稳定的单质。

标准摩尔反应焓的计算

Gaussian程序提供频率计算关键字（Freq）进行体系的热化学分析，因此，为了计算反应的焓变，我们可以直接在结构优化的同时加上Freq关键字。

计算体系反应焓变的基本流程：

  1. 优化分子结构，同时进行频率计算
  2. 在优化的结构基础上，用更大的基组进行单点能计算
  3. 计算标准摩尔反应焓变

在这个教程中，我们以1,3-丁二烯作为例子，计算反式构象转化为顺式构象时的标准摩尔反应焓变。

优化反式和顺式1,3-丁二烯的结构

反式丁二烯的结构优化与频率计算输入文件(trans-opt.com)：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

# b3lyp/6-31G* opt freq 
 
trans-OPT
 
0 1
 C                  0.61361546    0.29986053    0.00803267
 H                  1.14677921   -0.62784439    0.00803267
 H                 -0.45638454    0.29986053    0.00803267
 C                  1.28888977    1.47483782    0.00803267
 H                  0.75572603    2.40254274    0.00803267
 C                  2.82888977    1.47483782    0.00803267
 H                  3.36205352    0.54713290    0.00803267
 C                  3.51113317    2.64578231    0.00803267
 H                  2.98349070    3.57663859    0.00803267
 H                  4.58111429    2.63942501    0.00803267
 
!空白行不可忽略

计算:

1

$g16root/g16/g16 trains-opt.com

计算完毕，得结果文件trans-opt.log。

顺式丁二烯的结构优化与频率计算输入文件(cis-opt.com)：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

# b3lyp/6-31G* opt freq 
 
cis-OPT
 
0 1
 C                  0.61361546    0.29986053    0.00803267
 H                  1.14677921   -0.62784439    0.00803267
 H                 -0.45638454    0.29986053    0.00803267
 C                  1.28888977    1.47483782    0.00803267
 H                  0.75572603    2.40254274    0.00803267
 C                  2.82888977    1.47483782    0.00803267
 H                  3.36205352    2.18304610    0.60726338
 C                  3.51113317    0.58094088   -0.74831333
 H                  2.98349070   -0.12967314   -1.34957959
 H                  4.58111429    0.58579404   -0.74420697
 
!空白行不可忽略

计算:

1

$g16root/g16/g16 cis-opt.com

计算完毕，得结果文件cis-opt.log。

单点能计算

使用第一步优化得到的结构作为初始结构。在更大的基组下计算两个构象的单点能。

反式丁二烯的但电能计算输入文件(trans-sp.com)：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

# b3lyp/6-311+g(d,p)
 
trans-SP
 
0 1
 C                  1.10712900   -1.48465400    0.00000000
 H                  2.09998900   -1.03986800    0.00000000
 H                  1.05995700   -2.56940600    0.00000000
 C                 -0.00001100   -0.72877000    0.00000000
 H                 -0.97811200   -1.21104800    0.00000000
 C                  0.00000000    0.72875800    0.00000000
 H                  0.97811000    1.21101800    0.00000000
 C                 -1.10712300    1.48466300    0.00000000
 H                 -2.09999700    1.03991100    0.00000000
 H                 -1.05991600    2.56941400    0.00000000
 
!空白行不可忽略

计算:

1

$g16root/g16/g16 trains-sp.com

计算完毕，得结果文件trans-sp.log。

顺式丁二烯的但电能计算输入文件(cis-opt.com)：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17

# b3lyp/6-311+g(d,p)
 
cis-SP
 
0 1
 C                  1.54361300   -0.49345700   -0.07503600
 H                  1.17452700   -1.45548800   -0.42194500
 H                  2.61090500   -0.41878400    0.11204100
 C                  0.72819400    0.55430000    0.09993500
 H                  1.16348700    1.50921600    0.39654900
 C                 -0.72820400    0.55430800   -0.09993400
 H                 -1.16348400    1.50922200   -0.39655400
 C                 -1.54360900   -0.49346200    0.07503600
 H                 -1.17449600   -1.45548400    0.42193700
 H                 -2.61090600   -0.41881500   -0.11203000
 
!空白行不可忽略

计算:

1

$g16root/g16/g16 cis-sp.com

计算完毕，得结果文件cis-sp.log。

计算构象转化的标准摩尔反应焓

  1. 提取焓校正值

    2. 首先，分别打开第一步优化的两个构象的输出文件(trans-opt.log, cis-opt.log)，找到下面的输出部分（以反式构象为例），提取焓校正值。

    1
2
3
4
5
6
7
8

     Zero-point correction=                           0.085491 (Hartree/Particle)
 Thermal correction to Energy=                    0.090134
 Thermal correction to Enthalpy=                  0.091079
 Thermal correction to Gibbs Free Energy=         0.059049
 Sum of electronic and zero-point Energies=           -155.906647
 Sum of electronic and thermal Energies=              -155.902004
 Sum of electronic and thermal Enthalpies=            -155.901060
 Sum of electronic and thermal Free Energies=         -155.933089

    其中第3行，Thermal correction to Enthalpy,就是我们需要的焓校正值,后面计算要用到该值。

    其中第7行是我们需要计算的化合物的焓，但是我们用了一个更大的基组来计算电子能，理论上精度将更高，所以这里没有直接用该行给出的值。更多的该部分输出结果的解释请参考Joseph W. Ochterski博士的Thermochemistry in Gaussian一文[2]。

  2. 提取电子能量

    3. 再分别打开第二步两个构象单点能计算的输出文件（trans-sp.log, cis-sp.log），找到下面的输出部分（同样以反式构象为例），提取电子能。

    1
2
3
4

     Requested convergence on             energy=1.00D-06.
 No special actions if energy rises.
 SCF Done:  E(RB3LYP) =  -156.040758479     A.U. after   10 cycles
            NFock= 10  Conv=0.72D-08     -V/T= 2.0052

    其中第3行，以SCF Done开头为特征，E(RB3LYP)的值为电子能量。

  3. 计算标准摩尔反应焓变

    4. 利用同样的方法，我们提取了顺式和反式构象的焓校正值和电子能（提取的能量单位均为Hartree），并列于下表1：

    表1. 焓校正值、电子能与标准摩尔反应焓

    焓校正值 电子能 标准态下焓值 标准摩尔反应焓
    Trans-1,3-diene 0.091079 -156.040758479 -155.9496795 14.5 KJ/mol
    Cis-1,3-diene 0.090997 -156.035149002 -155.944152


    其中，

    标准态下焓值 = 焓校正值 + 电子能

    标准摩尔反应焓 = 产物的标准态焓值 – 反应物的标准态焓值

    1 Hartree = 2625.5 KJ/mol

    最终，我们得到反式构象转化为顺式构象的标准摩尔反应焓为14.5 KJ/mol。

标准摩尔生成焓的计算

接下来，我们以HCl分子为例，介绍如何计算分子的标准摩尔生成焓，上面已经提到标准摩尔生成焓是指由最稳定的单质生成该物质的焓变，也就是说，本质上它是一种特殊的标准摩尔反应焓，其在标准条件下的反应方程式如下：

1/2 H2 (g) + 1/2 Cl2 (g) → HCl (g)

在这个例子中，其最终的标准摩尔生成焓计算公式如下：

ΔfHmΘ = HmΘ(HCl) – 1/2HmΘ(H2) – 1/2HmΘ(Cl2)

在这里，我们将使用同样的方法基组和计算步骤计算三种分子的标准焓值，感兴趣的读者可以按照流程尝试练习，这里，我们将直接在表2中列出计算结果。

表2. 标准摩尔生成焓的计算

焓校正值（Hcor） 电子能(E) HmΘ ΔfHmΘ
H2 0.013450 -1.17956974143 -1.1661197 -100.8 KJ/mol
Cl2 0.004695 -920.409553181 -920.4048582
HCl 0.009980 -460.833843451 -460.8238635


将表2中的值代入上述公式，即可得到HCl分子的标准摩尔生成焓为 -100.8 KJ/mol。

原子化法计算标准摩尔生成焓

最后我们将介绍如何利用原子化法进行标准摩尔生成焓的计算，首先，假如分子M的化学式为AxBy，那么，生成分子AxBy的过程可以写成如下形式，也就是所谓的原子化形式：

AxBy → xA + yB

其中A和B是组成分子M的原子，x和y分别是组成分子M的A和B原子的个数。下面我们直接给出标准摩尔生成焓的计算公式：

ΔfHmΘ(M,298K) = HmΘ(M,298K) + xΔfHmΘ(A,0K) + yΔfHmΘ(B,0K) -xE(A,0K) -yE(B,0K)
–( x(HmΘ(A,298K) – HmΘ(A,0K)) + y(HmΘ(B,298K) – HmΘ(B,0K)))

其中，ΔfHmΘ(M,298K)是分子M的标准摩尔生成焓；HmΘ(T)是分子或原子在温度T时的焓；ΔfHmΘ(0K)是原子A或B在0K时的摩尔生成焓；E(0K)是原子A或B在0K时的电子能；需要指出的是，上面公式中红色字体高亮部分将采用实验测量值(表3)。

表3. 常见原子的自旋多重度、0K时的生成焓以及0K至298K温度校正（单位为kcal/mol）[2,3]

Elements 自旋多重度 ΔfHΘ(0K) HΘ(298 K) – HΘ(0 K)
H 2 51.63±0.001 1.01
Li 2 37.69±0.2 1.10
Be 1 76.48±1.2 0.46
B 2 136.2±0.2 0.29
C 3 169.98±0.1 0.25
N 4 112.53±0.02 1.04
O 3 58.99±0.02 1.04
F 2 18.47±0.07 1.05
Na 2 25.69±0.17 1.54
Mg 1 34.87±0.2 1.19
Al 2 78.23±1.0 1.08
Si 3 106.6±1.9 0.76
P 4 75.42±0.2 1.28
S 3 65.66±0.06 1.05
Cl 2 28.59±0.001 1.10

我们以呋喃为例，采用与上面相同的方法计算分子的焓值和原子的电子能，呋喃分子的原子化形式为：

C4H4O → 4C + 4H + O

将计算所得结果代入公式（由于呋喃由三种元素组成，相应的改变公式形式）：

ΔfHmΘ(M,298K) = HmΘ(furan,298K) + 4ΔfHmΘ(C,0K) + 4ΔfHmΘ(H,0K) + ΔfHmΘ(O,0K) -4E(C,0K) -4E(H,0K) -E(O,0K) – ( 4(HmΘ(C,298K) – HmΘ(C,0K)) + 4(HmΘ(H,298K) – HmΘ(H,0K)) + (HmΘ(O,298K) – HmΘ(O,0K)))

= -3899.862084 + 3958.326324 – 25.45574

= -30.0 KJ/mol

最终计算结果与其实验测量值-34.7 kJ/mol相差4.7 KJ/mol，通过提升计算精度，我们也能够更加逼近真实值。

参考文献

  1. Exploring Chemistry with Electronic Structure Methods, 3rd ed., Gaussian, Inc.: Wallingford, CT, 2015. J. B. Foresman and Æ Frisch
  2. Joseph W. Ochterski(2000). Thermochemistry in Gaussian. http://gaussian.com/thermo
  3. Curtiss, L. A.; Raghavachari, K.; Redfern, P. C.; Pople, J. A. Assessment of Gaussian-2 and Density Functional Theories for the Computation of Enthalpies of Formation. J. Chem. Phys. 1997, 106 (3), 1063–1079. https://doi.org/10.1063/1.473182.

采购Gaussian，请联系我们