第一性原理方法研究小芳香烴分子的基態性質

王金華 李澤朋

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/