第一性原理方法研究小芳香烃分子的基态性质

王金华 李泽朋

1.天津职业技术师范大学理学院，天津 300222

2.中国民航大学理学院，天津 300300

1 概述

芳香烃是含苯环结构的碳氢化合物，来源于煤、石油等，是化学工业中的重要原料[1]。碳黑、煤烟是有机燃料燃烧的主要产物。有机燃料的不完全燃烧不仅会产生常规污染物，如SO2、NO、CO，而且会产生致癌的污染物，像多环芳香烃（PAH），造成燃料浪费和环境污染。因此改善碳燃料燃烧的效率，减少其副产品的产量，成为实验和理论研究的热点之一[2，3]。苯是多环碳氢化合物的基本单元，多个苯环的结合可得到单层或者多层PAH化合物。所以可以通过研究小芳香烃及其自由基的性质，明确燃烧过程中的中间产物的性质，进而研究PAH分子的氧化反应来理解燃烧反应的潜在机制。

目前，有机燃料的燃烧机理已被广泛研究[1-3]。这些化合物的氧化是破坏它们的主要途径，通常氧化反应是定域的苯氧基参与反应，因此我们可以通过脱去氢、引入氧实现氧化，来研究发生氧化反应苯环及周围几个苯环的性质。理论方法可以有效地避免实验中反应条件复杂，难以得到有效数据的缺点，所以理论上可以有效地控制反应条件，选取合适的计算方法和参数对有机燃料的燃烧过程进行模拟。基于密度泛函（DFT）和二阶多体微扰理论（MP2）的第一性原理方法，能计算包含几十个苯环的PAH体系[4，5]。理论上对芳香烃性质的研究主要用基于密度泛函的杂化方法（BLYP、B3LYP、PW91、BMK、M06L等）、HF、MP2、CCSD方法。其中B3LYP方法计算得到的C-C键长、振动频率比HF方法得到的结果与实验结果符合得好[6]。HF近似方法得到的C-C键长比实验值小，振动频率比实验值高大约10%[7]。MP2微扰方法存在收敛性不一致的问题，但对芳香烃的研究能得到精确的结果[8]，但由于MP2方法计算花费大，计算小芳香烃分子时被广泛使用。

为了研究多环芳香烃的氧化反应机制，考虑到计算花费和效率，本文中应用基于DFT和MP2的第一性原理方法，研究小芳香烃（苯、萘）及其自由基的基态性质，如几何结构、振动频率和能量，并研究了氧原子在不同位置的五苯环（C19H10O）芳香烃的稳定性，为大PAH体系的研究及理解其燃烧反应的机制提供重要依据。

2 计算方法

以往的研究表明DFT和MP2方法能得到精确的分子结构[6，8]，本文计算中使用的方法包括：基于密度泛函 的 B3LYP/6-31G（d）、B3LYP/6-311++G（d，p）、B3LYP/6-311G（3df，3pd）、B3LYP/cc-pVDZ和B3LYP/cc-pVTZ方 法，MP2近 似 的RIMP2/ccpVDZ、RIMP2/cc-pVTZ 和 SOSMP2/cc-pVDZ、SOSMP2/cc-pVTZ方法；基于密度泛函的EDF1、EDF2函数使用文献中推荐的EDF1/6-31+G（d）、EDF2/ccpVTZ[9]。通过振动频率的计算，这些结构不存在虚频，优化能得到局域能量最小的结构。使用PW91、M06L、BMK方法对苯分子进行几何优化和振动频率的计算，优化后的结构与实验值有很大差异，振动频率结果显示存在虚频，所以不使用这些方法研究小芳香烃及其自由基的性质。以上计算都使用Q-Chem软件包中的相应程序完成[10]。苯的氧化反应路径为

多环芳香烃（PAH）的氧化反应与苯的氧化反应类似。本文中研究了苯和萘氧化反应中的基团以及小芳香烃，如苯酚、甲苯、萘酚和自由基，如甲基和羟基等。

3 结果和讨论

使用上述的DFT和MP2方法对所研究的结构进行几何优化，部分结构C-C键、C-H键、C-O键、O-H键键长的理论计算值和实验值在表1中列出，苯相应的键长、键角的平均绝对偏差如图1所示，实验值均来自计算化学对比与基准数据库（CCCBDB）[11]。从表1和图1中的数据可以看出，随着基底的增大，键长和键角与实验值的平均绝对误差呈无规律的变化，不能得到基底的大小对几何结构影响的趋势。DFT和MP2方法得到的C-C、C-H键与实验值的最大误差大约为1%，CCC角、HCC角与与实验值的最大误差为0.02%，都与实验值符合得很好。如苯分子，B3LYP/6-31G（d）方法得到的C-C键长与实验值（0.1397nm）最接近，相差0.00004nm。EDF2/cc-pVTZ方法得到苯的C-C键长为0.1386nm，与实验值相差最大为0.0011nm。B3LYP/6-311++G（d，p）方法得到的C-H键长的值与实验值（0.1084nm）最接近，相差0.00005nm。小芳香烃及其自由基的C-H键比C-C键弱，所以C-C键对几何结构的影响比C-H键强，可以根据C-C键的结果来评估各种方法对几何结构的影响。B3LYP/6-31G（d）和EDF2/cc-pVTZ分别给出与实验结果相比最优和最差的苯分子结构。实验上得到苯酚的C-C键长为0.1398nm，B3LYP/6-31G（d）方法得到的键长C-C键长（0.1397nm）与实验值最接近，相差0.0001nm。与苯类似，EDF2/ccpVTZ方法得到苯酚的C-C键长（0.1387nm）与实验值相差最大，为0.0011nm。B3LYP/6-31G（d）和EDF2/ccpVTZ方法计算得到的C-C键长与实验值的平均绝对偏差分别为0.1%和0.8%。如萘分子C-C键计算结果与实验结果相比较，B3LYP/cc-pVTZ和SOSMP2/cc-pVDZ方法分别给出最优和最差的结果。甲苯分子苯环上C-C键计算结果与实验结果相比较，B3LYP/6-311++G（d，p）和SOSMP2/cc-pVDZ方法分别给出最优和最差的结构。对所研究的小芳香烃及其自由基，有的自由基没有实验数据，假定给出自由分子最优和最差结构的方法，同样能给出自由基的最优和最差结构，比如，B3LYP/6-31G（d）和EDF2/cc-pVTZ方法给出苯基、苯氧基、环戊二烯基最优和最差的结构。研究结果表明，对所研究的结构，C-C、C-H键、键角的误差很小，所使用的DFT和MP2方法都可以很精确地研究这些体系的结构。DFT中的B3LYP方法得到的结果与实验结果符合得最好。虽然还有更精确的方法，像量子蒙特卡罗方法

（QMC）、耦合团簇法（CC），但考虑到研究大PAH体系时的计算效率，B3LYP方法是在一定的计算精度和效率内的最佳选择之一。

表2 CCSD（T）/6-31G（d）方法计算得到的最优与最差结构之间的相对能Er其中茚基使用CCSD/6-31G（d）方法

表1 DFT方法和MP2方法计算得到的小芳香烃及自由基的键长（10-1nm）及实验值。

图1 DFT和MP2方法计算得到苯的（a）C-C、（b）C-H、（c）CCC角、（d）HCC角与实验值的平均绝对偏差。基底从左到右分别为：6-31G（d）、6-31+G（d）、6-311++G（d，p）、6-311G（3df，3pd）、cc-pVDZ 和 cc-pVTZ

小芳香烃分子及其自由基的几何优化得到的C-C和C-H键的误差都很小，然而C-H键比C-C键弱，所以C-C键的变化对能量的影响较大。计算得到的与实验结果比较最优和最差的几何结构，使用比B3LYP更精确的耦合团簇方法CCSD（T）（除了茚基使用CCSD）来计算它们的能量，以得到几何结构的不同对能量的影响。CCSD（T）使用6-31G（d）和cc-pVTZ基底计算得到的能量差别很小，在我们的计算中使用小基底6-31G（d），相应的结果在表2中列出。计算得到，甲基、苯甲基和甲苯最优和最差的几何结构之间能量差别很小，分别为0.03kJ/mol、0.42kJ/mol和0.51kJ/mol；萘酚和萘氧基，最优和最差的几何结构之间能量差别相对较大，分别是6.62和8.05kJ/mol。对萘和羟基，能量不同分别是-1.90kJ/mol和-0.74kJ/mol，最优的几何结构比最差的结构的能量高。CCSD（T）方法没有得到茚基的基态能量，计算中使用CCSD/6-31G（d）得到相对能量为1.48kJ/mol。所有小芳香烃及其自由基，DFT和MP2方法优化得到的最优和最差结构间相对能的最大差异大约是8kJ/mol，基底越大，结果越精确，但是对分子能量的影响不是很大，因此可以用B3LYP方法优化后得到的最优结构来研究小芳香烃分子及自由基的基态能量，同时也可以研究它们的其他性质，像热力学特性。B3LYP方法在计算精度和效率方面有优势，是模拟多环芳香烃（PAH）化学性质的最佳方法之一，可以推广到大PAH的研究中去。

使用B3LYP方法研究五苯环（C19H10O）中氧原子的位置对系统稳定性的影响，B3LYP/6-311G（d，p）方法几何优化的结构如图2所示。表3列出了B3LYP使用6-31G（d）、6-311G（d，p）、cc-pVTZ基底计算得到的不同结构的相对能Er。基底越大，能量越低，但三种不同的基底对能量的影响不大。比较（a）-（f）结构的能量得到，氧原子在（d）位置时的能量最低，结构最稳定，相应的C=O键最短为0.1224nm。氧原子在（b）、（e）位置时的能量较高，结构不稳定，相应的C=O键为0.1240nm。C=O键越短，体系越稳定。六种结构相对稳定性为：（d）＞（a）＞（c）＞（f）＞（e）＞（b）。从以上分析可以看出，B3LYP方法适用于研究五苯环体系的性质，进一步可以研究大的多环芳香烃甚至石墨烯体系的性质。

图2 氧原子在不同位置的五苯环结构（a）、（b）、（c）、（d）、（e）、（f）

表3 B3LYP/6-311G（d，p）优化后的五苯环结构，B3LYP使用6-31G（d）、6-311G（d，p）、cc-pVTZ基底计算得到的相对能Er。

4 结束语

采用DFT和MP2方法研究小芳香烃分子及其自由基的几何结构和能量。研究表明，B3LYP方法能给出与实验结果非常符合的基态结构。同时研究了氧原子在不同位置的五苯环体系稳定性，得到氧在（d）位置时体系最稳定。B3LYP方法是模拟小芳香烃化学性质的最佳方法之一，可以推广到大PAH的研究中去。

[1]RICHTER. H，HOWARD J B. Formation of polycyclic aromatic hydrocarbons and their growth to soot—a review of chemical reaction pathways [J]. Prog. Energ. Combust. Sci.2000，26（4-6）: 565-608.

[2]ABDULLAH H H，KUBBA R M，SHANSHAL M.Vibration Frequencies Shifts of Naphthalene and Anthracene as Caused by Different Molecular Charges [J]. Z. Naturforsch.2003，58a，645-655.

[3]FIROUZI R，ZAHEDI M. Polyacenes electronic properties and their dependence on molecular size [J]. J. Mol. Struct. :THEOCHEM 2008，862（1-3）: 7-15.

[4]MIANI A，CANE E，PALMIERI P，et al. Experimental and theoretical anharmonicity for benzene using density functional theory [J]. J. Chem. Phys. 2000，112（1）: 248-259.

[5]BLANKSBY S J，ELLISON G B. Bond Dissociation Energies of Organic Molecules，Acc. Chem. Res. 2003，36（4）: 255-263.

[6]HAJGATO B，DELEUZE M S，TOZER J，et al. A benchmark theoretical study of the electron affinities of benzene and linear acenes [J]. J. Chem. Phys. 2008，129，084308.

[7]RIVES J T，JORGENSEN W L. Performance of B3LYP Density Functional Methods for a Large Set of Organic Molecules [J]. J. Chem. Theor. Comput. 2008，4（2）:297-306.

[8]ASHVAR C S，DEVLIN F J，STEPHENS P J. Molecular Structure in Solution: An ab Initio Vibrational Spectroscopy Study of Phenyloxirane [J]. J. Am. Chem. Soc. 1999，121（12），2836-2849.

[9]MERRICK J P，MORAN D，RADOM L. An Evaluation of Harmonic Vibrational Frequency Scale Factors [J]. J. Phys.Chem. A，2007，111（45）: 11683-11700.

[10]SHAO，Y，MOLNAR L F，JUNG Y，et al. Advances in quantum chemical methods and algorithms in the Q-Chem 3.0 program package [J]. Physical Chemistry Chemical Physics 2006，8（27）: 3172-3191.

[11]R. D. Johnson，NIST Computational Chemistry Comparison and Benchmark Database： NIST Standard Reference Database，Number 101 [DB/OL]. （2013-8）[2005-5]. http://cccbdb. nist. gov/