劉 珉,魏賢鳳
高能炸藥CL-20分子結構的理論模擬方法探究
劉 珉,魏賢鳳
(西南科技大學國防科技學院,四川 綿陽,621010)
采用量子化學方法,在HF/6-31G, DFT-B3LYP/6-31G 和MP2/6-31G 基組水平下全優(yōu)化ε-CL-20分子的結構,對ε-CL-20的結構(包括鍵長、鍵角、二面角)分析得出B3LYP/6-31G水平下的理論值更接近實驗值。采用B3LYP 方法,不同基組3-21G,6-31G, 6-31++G, 6-311++G,6-31++G**水平下全優(yōu)化ε-CL-20分子,基組選取基本不影響計算結果, 表明B3LYP/6-31G 基組水平下的計算結果能夠滿足ε-CL-20的結構優(yōu)化要求。采用HF、B3LYP、MP2方法不同基組計算ε-CL-20原子的凈電荷分布, 結果表明充分考慮電子相關性的B3LYP 方法計算結果最合理。
炸藥;CL-20;量子化學;B3LYP/6-31G;分子結構;HF;MP2
六硝基六氮雜異伍茲烷(HNIW,俗稱CL-20)由美國的 Nielsen等[1-2]在 1987年首次合成,是一種最有發(fā)展?jié)摿Φ母吣苊芏然衔镏?。CL-20在常溫下有4種晶型,即α、β、γ和ε,其中ε-CL-20[3]的分子對稱性最高,熱力學性能最穩(wěn)定。量子化學[4]計算方法在預測材料性能方面已得到很大的應用,Hartree-Fock從頭計算方法、密度泛函理論(DFT)已被應用于計算晶體的結構和性質[5-7]。肖鶴鳴等人[8]首次運用DFT方法研究CL-20 4種晶型的能帶和電子結構,預測感度大小,驗證得知ε-CL-20感度最低。
本文運用Hartree-Fock從頭計算方法、DFT、MP2方法在6-31G基組水平上對CL-20結構做全優(yōu)化計算,并在B3LYP方法下采用不同的基組水平對ε-CL-20分子做全優(yōu)化計算,分析不同計算水平下ε-CL-20的結構和性能,為今后用量子化學方法研究CL-20提供理論支持。
CL-20常溫下有4種晶型,即α、β、γ和ε。主要區(qū)別在于它們的硝基在5個氮環(huán)平面的取向(即內(nèi)環(huán)和外環(huán)之間的差異)上是不同的。其中,ε型熱力學性能最穩(wěn)定,晶型如圖1所示。本文采用CL-20最佳幾何結構[9],用Gaussian View搭建模型Gaussian09程序全優(yōu)化計算,采用HF、MP2、B3LYP方法在6-31G水平上全優(yōu)化分子結構,優(yōu)化結果收斂正常且無虛頻,證明結構是最穩(wěn)定狀態(tài);并采用B3LYP在3-21G、6-31G、6-31++G、6-311++G、6-31++G**水平上全優(yōu)化計算分子結構,分析優(yōu)化后的結構和性能。
圖1 ε-CL-20的晶型
圖2為ε-CL-20分子的原子編號和優(yōu)化后的幾何示意圖。由圖2可見,CL-20是一多氮多環(huán)籠形硝胺分子,其基本結構為一籠形的剛性異伍茲烷,6個硝基分別連接于6個橋氮原子上。
圖2 ε-CL-20的結構圖
表1~3 分別列出ε-CL-20的全優(yōu)化幾何構型參數(shù)鍵長、鍵角、二面角,由于此分子結構具有一定的對稱性,故表中只列出部分數(shù)據(jù)。由表1可以看出,HF方法計算的鍵長比實驗值偏短,平均相差0.030nm。MP2方法計算的鍵長與實驗值接近,平均相差0.014nm。B3LYP方法計算的鍵長也比實驗值偏短,平均相差0.001nm。
表1 不同方法在6-31G基組下計算ε-CL-20的鍵長
Tab.1 The calculated bond lengths of ε-CL-20 at the level of HF/ 6-31G, MP2/ 6-31G and B3LYP/ 6-31G
注:誤差值=計算值-參考值
表2 不同方法在6-31G基組下計算ε-CL-20的鍵角
Tab.2 The calculated bond angles of ε-CL-20 at the level of HF/ 6-31G, MP2/ 6-31G and B3LYP/ 6-31G
根據(jù)表1中的數(shù)據(jù),3種方法計算的鍵長相差不大,與實驗值非常接近。在小基組下,HF只是考慮了電子自旋,而沒有考慮電子相關能,所以鍵長比較短。但是當基組比較小時,HF的鍵長計算值反而比B3LYP和MP2都精確,這是因為小基組的函數(shù)不能滿足B3LYP和MP2法的計算要求,所以出現(xiàn)比較大的誤差。而且在任何計算方法下,C-H鍵的鍵長誤差是最小的,而N-N鍵和N-O鍵的誤差比較大,說明在CL-20分子中,C-H上電子分布比N-N鍵和N-O鍵上的電子分布均勻。另外由于C-N單雙鍵長分別是1.471nm和1.273nm,所以C-N鍵為單鍵,所有C原子為sp3雜化。而N-N單鍵長度為1.47nm,N-N雙鍵的長度為1.25nm,所以N-N介于單雙鍵之間。
由表2可以看出,HF方法計算的鍵角比實驗值偏大,平均相差0.659°。MP2法計算的鍵角與實驗值接近,平均相差0.221°。B3LYP方法計算的鍵角也比實驗值偏大,平均相差0.065°。由表2中數(shù)據(jù)可知,使用3種計算方法對CL-20鍵角的影響不一樣,只有B3LYP方法最接近實驗值。其中,以C為中心原子的鍵角計算值接近于110°,而以N為中心原子的鍵角計算值接近于119°,由于ε-CL-20分子的中心對稱性,可以知道環(huán)上的所有C原子采取sp3雜化,環(huán)上N原子和硝基N原子均采取sp2雜化。
表3 不同方法在6-31G基組下計算ε-CL-20的二面角
Tab.3 The calculated dihedral angles of ε-CL-20 at the level of HF/ 6-31G, MP2/ 6-31G and B3LYP/ 6-31G
由表3可以看出,HF方法計算的二面角比實驗值偏小,平均相差2.4°。MP2法計算的鍵長與實驗值接近,平均相差1.2°。B3LYP方法計算的鍵長也比實驗值偏小,平均相差0.60°。二面角的計算值有一定的偏差,但是MP2方法和B3LYP的計算值更接近,ε-CL-20的分子構型和文獻相符。
綜上所述,使用MP2,HF,B3LYP這3種方法對ε-CL-20結構進行優(yōu)化,計算得到的結果都是在實驗值誤差范圍之內(nèi)。但是,考慮到精度因素,B3LYP方法更接近于實驗值,使用B3LYP方法優(yōu)化ε-CL-20結構更合適。
基于2.1研究可知采用B3LYP 方法優(yōu)化出的分子幾何構型結果更好。因此,采用B3LYP 方法,取用不同的基組:3-21G,6-31G,6-31++G,6-311++G, 6-31++G**,對結構再次進行優(yōu)化計算。結果分別見表4~6。
表4 B3LYP方法不同基組下ε-CL-20的計算鍵長
Tab.4 The calculated bond lengths of ε-CL-20 at different levels with B3LYP method
表5 B3LYP方法不同基組下ε-CL-20的計算鍵角
Tab.5 The calculated bond angle of ε-CL-20 at different levels with B3LYP method
由表4可以看出,在B3LYP的方法下,除3-21G和6-311++G**基組的計算值以外,其他基組的鍵長計算值相差都不大。在6-31++G基組下,N-O鍵的鍵長與實驗值平均相差0.005nm,其他鍵長基本沒有出入。除基組 3-21G和 6-31++G**,其他幾個基組的變化不會影響計算的結果,接近于實驗值。由表5可以看出,基組的變化對鍵角的影響不大,以N為中心原子的鍵角約為120°,以C為中心原子的鍵角約為110°,都在正常范圍內(nèi),接近實驗值。從表6可以看出,基組的變化對二面角的影響不大,計算出的二面角在正常范圍內(nèi),因此結構基本不變。
綜上所述,相同方法、不同基組計算ε-CL-20分子時結果并沒有特別大的差異。但是,與實驗值最接近的計算結果是6-311++G基組水平。因為晶體中分子受晶格能束縛和周圍分子的影響,而且計算的結果為理想氣體狀態(tài)下的分子結構,所選用的基組也有一定近似性,所以,優(yōu)化計算所得構型與晶體狀態(tài)下的分子結構不可能完全一致,計算值與晶體數(shù)據(jù)之間是存在一定偏差的。因此綜合考慮,選取B3LYP/6-31G優(yōu)化的結構基本可以滿足需求。
表6 B3LYP方法不同基組下ε-CL-20的計算二面角
Tab.6 The calculated dihedral angle of ε-CL-20 at different levels with B3LYP method
不同方法全優(yōu)化-CL-20分子的時間圖見圖3。
圖3 不同方法優(yōu)化ε-CL-20分子的時間柱形圖
由圖3可以看出,在相同基組、不同方法下,計算的時間長短順序為HF 不同方法不同基組下CL-20的單點能計算如表7所示。 表7 不同方法計算CL-20的單點能(hartree) Tab.7 The calculated single point energy of ε-CL-20 with different methods 由表7可以看出,不同方法不同基組下CL-20的能量相差不大,與實驗值一致。HF法因為沒有考慮電子相關,MP2法只考慮了60%的相關能,整個計算結果的能量都偏低[10]。在B3LYP方法中3-21G基組與其他基組計算結果相差較大,而其他基組呈現(xiàn)出基組越大,能量越大的趨勢。這是因為大基組加入了極化、彌散等條件,分子中更多的能量被考慮到,而且大基組分為多個基函數(shù),使誤差降到更低[11]。綜上所述,小基組的計算值偏小,大基組考慮了更多的相關能后,能量增加,總體來說CL-20的能量變化不大,都在誤差允許范圍內(nèi)。通過比較HF、MP2、B3LYP與實驗值,發(fā)現(xiàn)B3LYP的計算結果最接近實驗值。而采用B3LYP方法不同基組3-21G、6-31G、 6-31++G、6-311++G、6-31++G**水平下,基組越大,計算精度越高,但時間也越長。綜合考慮B3LYP/ 6-31G為ε-CL-20結構計算的最優(yōu)方法。 電荷集居數(shù)分布是由電子波函數(shù)決定的分子中的電子在各原子上和原子間的分布密度。通過不同的方法、不同的基組計算了ε-CL-20上各原子電荷的分布情況。由于電荷分布具有中心對稱性,表8中只列出了部分原子的電荷分布。 由表8可以看出,通過不同的方法、不同基組計算ε-CL-20上各原子電荷的分布情況差別較大。在同一基組(6-31G)、不同方法水平下,亞硝基上的N原子帶正電荷,O原子帶負電荷且相差不大,H原子帶正電荷,環(huán)上的C原子帶正電荷,環(huán)上的N原子帶負電荷。 采用B3LYP方法計算不同基組3-21G、 6-31G、 6-31++G、 6-311++G、6-31++G**水平下的原子凈電荷。除3-21G和 6-31G基組外,其它基組的計算結果是環(huán)上的C原子帶負電荷,環(huán)上的N原子為正電荷,這與N原子上的分析結果和NO2有很大的不同。電荷不規(guī)則,O原子電荷接近0。因此,采用B3LYP方法在 3-21G、 6-31++G、 6-311++G和 6-31++G**基組水平下均不適合計算凈電荷的分布?;M3-21G太小,會產(chǎn)生較大的誤差;而基組 6-31++G、 6-311++G、6-31++G**可能是因為增加了極化和彌散函數(shù)使結果異常;基組6-31G的計算結果最合理。 總之,采用不同方法相同基組與參考值的比較,得出不同方法的精度順序從低到高依次為HF、MP2、B3LYP;而在B3LYP方法不同基組水平下分析計算結果得出基組6-31G計算準確。 表8 不同水平下ε-CL-20原子電荷的分布 Tab.8 The net charges distribution of ε-CL-20 atoms at different levels (1)對于高能炸藥ε-CL-20,充分考慮電子相關性的B3LYP方法計算的結果比HF方法計算的結果要更準確,并且耗費的時間比MP2要少。(2)綜合考慮時間和計算精度兩方面因素,用B3LYP/6-31G 水平基本上就可以滿足對ε-CL-20 結構優(yōu)化和單點能計算的要求。 [1] Nielsen A T.Cagedpolynitramine compound:US,5693794 [P]. 1997-12-02. [2] 李文祥,趙省向,邢曉玲.CL-20的合成、性能和應用研究進展[J].山西化工, 2015, 35(2):28-30. [3] 文國. CL-20晶型研究及CL-20/TATB共晶的分子動力學模擬[D].太原:中北大學, 2014. [4] 龍威.理論研究中的量子化學計算方法[J].寧夏師范學院學報, 2010,31(3):43-47. [5] Zhu W H,Xiao J J,Xiao H M.Comparative first-principles study of structural and optical properties of alkali metal azides[J]. J Phys Chem B, 2006(110):9 856-9 862. [6] Zhu W H,Xiao J J,Xiao H M.Density functional theory study of the structural and optical propertirs of lithium azide[J].Chem Phys Lett,2006(422):117-121. [7] Zhu W H,Xiao H M.Ab. Initio study of energetic solids:cupric azide,mercuric azides,and lead azide[J].J Phys Chem B,2006(110):18 196-18 203. [8] Xu X J,Zhu W H,Xiao H M.DFT studies on the four polymorphs of crystalline ε-CL-20 and the influences of hydro- static pressure on CL-20 crystal[J].J Phys Chem B,2007(111):2 090-2 097. [9] 趙信歧,施倪承.ε-六硝基六氮雜異伍茲烷的晶體結構[J].科學通報,1995(40):2 158-2 160. [10] Zhu K,Achord P D, Zhang X, et al. Highly effective pincer-ligated iridium catalysts for alkane dehydrogenation DFT calculations of relevant thermodynamic, kinetic, and spectroscopic properties.[J]. Journal of the American Chemical Society, 2004, 126(40):13 044-53. [11] Ciofini I, Daul C A. DFT calculations of molecular magnetic properties of coordination compounds[J].Coordination Chemi- stry Reviews, 2003, 238(02):187-209. Study on Methods for Theoretical Modeling the Structure of High Explosive CL-20 LIU Min,WEI Xian-feng (Southwest University of Science and Technology, Mianyang, 621010) Quantum chemistry methods were used to optimize the structure of ε-CL-20 at the level of HF/6-31G, DFT-B3LYP/6-31G and MP2/6-31G respectively. The structure of ε-CL-20, such as bond length, bond angle, dihedral angle,was analyzed theoretically.The theoretical value at the level of B3LYP/6-31G accords with the experimental value. The effects of B3LYP/3-21G, B3LYP/6-31G, B3LYP/6-31++G, B3LYP/6-311++G and B3LYP/6-31++G** methods on the calculated results of ε-CL-20 are small, which shows that the calculated results at the level of B3LYP/6-31G can meet the optimization of ε-CL-20. The atom charges distribution of ε-CL-20 at different levels were calculated quantitatively. The results show that the B3LYP method in which electron correlation effect has been sufficiently considered is more reasonable. Explosives;CL-20;Quantum chemistry;B3LYP/6-31G;Molecular structure;HF;MP2 TQ564 A 10.3969/j.issn.1003-1480.2019.06.007 1003-1480(2019)06-0025-05 2019-09-17 劉珉(1986 -),女,碩士研究生,主要從事火工品設計、含能材料模擬研究。 中國國家自然科學基金會財政支持(項目批準立項號11602239)。2.4 CL-20的電荷集居數(shù)(Mulliken)的分析
3 結論