劉 珉，魏賢鳳

高能炸藥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提供理論支持。

1 計算方法

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 結果與討論

2.1 不同q量化方法優(yōu)化ε-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.2 不同基組水平下優(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

2.3 CL-20的計算時間和單點能分析

不同方法全優(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)方法。

2.4 CL-20的電荷集居數(shù)（Mulliken）的分析

電荷集居數(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

3 結論

（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)。