劉英哲, 來蔚鵬, 尉 濤, 葛忠學(xué), 徐 濤, 駱艷嬌, 尹世偉

(1. 西安近代化學(xué)研究所氟氮化工資源高效開發(fā)與利用國家重點實驗室, 陜西 西安 710065;2. 陜西師范大學(xué), 陜西 西安 710062)

1 引 言

全氮材料具有能量高、無污染的優(yōu)點,是潛在的新型高能量密度材料[1-5]。然而,全氮材料的合成尚處于實驗室探索階段,難以獲得其實際樣品,如何獲取全氮材料在世界范圍內(nèi)仍然是一個巨大的挑戰(zhàn)。因此,對于全氮材料的能量性質(zhì),特別是爆轟性能,只能通過理論方法進(jìn)行預(yù)測。

晶體密度是決定全氮材料爆轟性能的重要參數(shù),其預(yù)測方法主要有兩類: 一是基于分子體積的半經(jīng)驗方法[6-10]; 一是基于晶體體積的分子力學(xué)方法[11-12]。前者缺乏實驗數(shù)據(jù)難以獲得準(zhǔn)確的經(jīng)驗參數(shù),且忽略了分子間相互作用對晶體堆積結(jié)構(gòu)的影響,在原理上不夠精確。后者缺乏完全適用于全氮結(jié)構(gòu)的力場參數(shù),無法準(zhǔn)確描述全氮分子間的相互作用,致使晶體結(jié)構(gòu)預(yù)測存在較大偏差。

目前,國際上關(guān)于全氮材料晶體密度預(yù)測研究中,最有代表性和被廣泛引用的報道來自于——英國QinetiQ公司[13]和瑞典國防研究院FOI[14]。二者均采用分子力學(xué)方法預(yù)測了5種籠形全氮結(jié)構(gòu)的晶體密度(見圖1),但預(yù)測結(jié)果卻存在較大差異。例如,對于N4(Td)結(jié)構(gòu),QinetiQ預(yù)測的晶體密度為1.752 g·cm-3,而FOI的預(yù)測結(jié)果為2.3 g·cm-3。

圖1典型籠形全氮分子結(jié)構(gòu)示意圖

Fig.1Schematic diagram of typical all-nitrogen molecular structures with cage type

為了準(zhǔn)確預(yù)測全氮材料的晶體密度,本研究基于量子化學(xué)計算方法,建立了適用于全氮結(jié)構(gòu)的力場參數(shù)。與傳統(tǒng)分子力場相比,該力場具有“特異性”,即每個全氮結(jié)構(gòu)均具有一套獨立的參數(shù),尤其是范德華參數(shù); 而在傳統(tǒng)分子力場中,通常是一類分子采用一套相同的參數(shù)。例如,對于籠形全氮結(jié)構(gòu)均采用相同的范德華參數(shù),這將難以區(qū)分全氮結(jié)構(gòu)差異對晶體密度的影響,勢必帶來較大誤差。

通過建立特異性力場參數(shù),本研究共預(yù)測了20種全氮材料的晶體結(jié)構(gòu),基于晶體體積計算了晶體密度。為了與國外文獻(xiàn)相比較,僅列出了5種籠形全氮結(jié)構(gòu)的晶體密度預(yù)測結(jié)果,重點介紹方法的建立。

2 計算方法

2.1 分子內(nèi)力場參數(shù)

采用密度泛函理論B3LYP/6-31G(d)方法[15-16],通過Gaussian程序[17]優(yōu)化全氮分子的幾何構(gòu)型,在最優(yōu)構(gòu)型的基礎(chǔ)上進(jìn)行振動頻率分析,借助VMD程序[18]中Paratool插件將Hessian矩陣中的特征振動模式分解投影到內(nèi)坐標(biāo)上,得到鍵長、鍵角、二面角的平衡值及相應(yīng)的力常數(shù)等分子內(nèi)力場參數(shù)。分子內(nèi)勢能函數(shù)形式如下:

Eintra=∑Ebond+∑Eangle+∑Edihedral

(1)

(2)

(3)

(4)

式中,Ebond,Eangle和Edihedral分別為化學(xué)鍵、鍵角與二面角勢能項;b0,θ0,φ0分別為鍵長、鍵角及二面角的平衡值;Kb,Kθ,Kφ為相應(yīng)的力常數(shù);n為多重度參數(shù)。建立分子內(nèi)力場參數(shù)的目的是為了保證全氮結(jié)構(gòu)的合理性,尤其是在分子間范德華參數(shù)優(yōu)化過程中,需要避免全氮結(jié)構(gòu)發(fā)生較大偏差。

2.2 分子間力場參數(shù)

分子間勢能分為靜電與范德華相互作用能,分別由庫侖定律和Lennard-Jones函數(shù)描述,具體函數(shù)形式如下:

Einter=∑Eelec+∑EvdW

(5)

(6)

(7)

式中,Eelec和EvdW分別為靜電和范德華相互作用勢能項;i和j為不同分子上的氮原子;qi和qj為原子點電荷;D0為介電常數(shù);rij為兩個原子之間的距離;ε為勢阱深度;σ是作用勢為0時原子間的距離。

全氮分子只有一種原子類型N,其電負(fù)性是一致的,對于對稱性高的結(jié)構(gòu)(如文中所列籠形結(jié)構(gòu)),原子點電荷為0。對于對稱性較差的全氮結(jié)構(gòu),由于氮原子周圍化學(xué)環(huán)境不盡相同,因此具有微弱的靜電相互作用。總之,全氮分子間作用主要由范德華參數(shù)主導(dǎo),建立準(zhǔn)確的范德華參數(shù)對于全氮晶體密度的預(yù)測至關(guān)重要。

首先,采用B3LYP/6-31G(d)方法擬合全氮分子靜電勢得到ESP電荷,將其指定為原子點電荷,由庫侖定律計算分子間靜電相互作用能。對于范德華參數(shù),采用文獻(xiàn)[19]的方法進(jìn)行優(yōu)化,主要原理為: 在盡可能窮盡全氮分子二聚體結(jié)構(gòu)的條件下,優(yōu)化范德華參數(shù)使分子力學(xué)方法計算的二聚體結(jié)合能與量子化學(xué)計算的二聚體結(jié)合能均方根誤差(RMSE)最小。RMSE定義如下:

(8)

為了避免二聚體抽樣的人為性與任意性,采用Monte-Carlo方法隨機(jī)產(chǎn)生一定密度下含有大量全氮分子的無定形結(jié)構(gòu),以任意兩個全氮分子之間的質(zhì)心距離為判據(jù)進(jìn)行二聚體的抽樣。無定形結(jié)構(gòu)由Material Studio軟件[20]中Amorphous Cell模塊構(gòu)建,密度設(shè)置為0.6～0.8 g·cm-3,范德華參數(shù)采用Dreiding力場[21]默認(rèn)參數(shù),其余力場參數(shù)為經(jīng)優(yōu)化后的新參數(shù)。

對每個抽取的二聚體結(jié)構(gòu),均采用MP2/cc-pvtz[22-23]方法計算結(jié)合能Ebind:

Ebind=EAB-EA-EB

(9)

式中,EAB,EA和EB分別為二聚體AB、二聚體中單體A與B的能量,J。在分別得到量子化學(xué)和分子力學(xué)計算的二聚體結(jié)合能之后,使用Quasi-Newton 算法對范德華參數(shù)進(jìn)行優(yōu)化。

2.3 晶體堆積結(jié)構(gòu)

基于建立的全氮分子特異性力場參數(shù),使用Materials Studio程序[20]中Polymorph模塊,將全氮分子按照晶體空間點群進(jìn)行堆積,預(yù)測全氮分子的晶體結(jié)構(gòu),取晶格能最低的晶型,根據(jù)晶胞參數(shù)計算晶體密度。

3 結(jié)果與討論

全氮分子的晶體密度主要取決于分子間相互作用,因此范德華參數(shù)的準(zhǔn)確程度對于晶體密度預(yù)測至關(guān)重要。由范德華參數(shù)的優(yōu)化過程可知,影響范德華參數(shù)的因素很多,諸如二聚體抽樣個數(shù),二聚體結(jié)合能是否校正,范德華參數(shù)種類等,為了得到相對穩(wěn)定的范德華參數(shù),分別對以上影響因素作了考察。另外,由于Dreiding力場中范德華參數(shù)也采用了Lennard-Jones(12-6)函數(shù)形式,且在二聚體構(gòu)建中作為初始值賦予全氮結(jié)構(gòu),因此,本研究優(yōu)化的范德華參數(shù)與Dreiding力場參數(shù)也進(jìn)行了比較。需要注意的是: Dreiding力場中N原子的范德華參數(shù)具有唯一數(shù)值,即ε=0.3238 kJ·mol-1、σ=3.6621 ?,為了表達(dá)簡潔,在余下的表格中不再提及,而僅列出本文優(yōu)化所得的范德華參數(shù)。

3.1 抽樣個數(shù)

以200個N4(Td)分子構(gòu)建的無定形結(jié)構(gòu)為例,密度分別設(shè)置為0.8和0.6 g·cm-3,抽取任意兩個質(zhì)心距離小于15 ?的二聚體結(jié)構(gòu),抽樣個數(shù)分別為4382和3231。表1為不同抽樣個數(shù)下優(yōu)化的N4(Td)范德華參數(shù)。從表1中可以看出,兩種抽樣個數(shù)下優(yōu)化的范德華參數(shù)相差不大,RMSE提高不到1%,說明基于無定形結(jié)構(gòu)的抽樣策略能夠得到比較穩(wěn)定的范德華參數(shù)。另外,與Dreiding力場參數(shù)相比,優(yōu)化后ε變大、σ變小,說明N4(Td)二聚體相互作用更強(qiáng)、結(jié)合更加緊密,且RMSE下降了約一個數(shù)量級,說明該參數(shù)能夠準(zhǔn)確地描述N4(Td)二聚體的分子間相互作用。

表1抽樣個數(shù)對范德華參數(shù)的影響

Table1Effect of sampling number on the van der Waals (vdW) parameters

samplenumbervdWparametersε/kJ·mol-1σ/?RMSEDreidingthiswork43820．67993．37620．2110．03432310．69123．36100．2320．035

Note： RMSE is rootmean square error.

3.2 結(jié)合能校正

通過量子化學(xué)計算全氮二聚體結(jié)合能時,由于存在基組重疊效應(yīng),采用式(9)得到的二聚體結(jié)合能一般要大于真實值,需要進(jìn)行校正。以N10(D2h)為例,采用counterpoise方法[24]對結(jié)合能進(jìn)行了校正,并考察了其對范德華參數(shù)的影響。如表2所示,結(jié)合能經(jīng)過校正后,優(yōu)化的范德華參數(shù)與Dreiding力場參數(shù)的RMSE相差不大,說明Dreiding力場能夠較好地描述環(huán)型全氮分子,但對籠形全氮分子的預(yù)測精度較差(見表1)。如果結(jié)合能不經(jīng)過校正,Dreiding力場參數(shù)的RMSE則從0.358增加到0.686。對比校正前后的范德華參數(shù)可以發(fā)現(xiàn),未校正的ε較大、σ較小,即分子間相互作用更強(qiáng),從而高估了二聚體結(jié)合能。因此,對結(jié)合能進(jìn)行校正是必要的。

表2結(jié)合能校正及參數(shù)個數(shù)對范德華參數(shù)的影響

Table2Effect of binding energy correction and parameter numbers on the van der Waals parameters

structureEbindcorrectionvdWparametersε1/kJ·mol-1σ1/?ε2/kJ·mol-1σ2/?RMSEDreidingthisworkYESNO0．36773．68990．3580．3320．49373．42400．35943．73140．3580．3310．58743．53450．6860．3460．46653．28740．67453．51910．6860．331

3.3 參數(shù)種類

基于全氮結(jié)構(gòu)對稱性,可以根據(jù)化學(xué)環(huán)境定義不同的氮原子類型。例如,對于N4(Td)分子,每個氮原子周圍化學(xué)環(huán)境相同,只需定義一種原子類型; 對于N10(D2h)分子,則可以進(jìn)一步分為兩種原子類型,即環(huán)上原子與橋原子。因此,以N10(D2h)為例,考察了原子類型種類對范德華參數(shù)的影響。如表2所示,無論結(jié)合能校正與否,雙原子類型與單原子類型相比,RMSE相差不大,準(zhǔn)確程度提高有限,說明采用單原子類型已經(jīng)能夠很好地描述全氮分子間相互作用。

綜上,采用無定形結(jié)構(gòu)抽樣方法可以得到相對穩(wěn)定的范德華參數(shù),二聚體結(jié)合能的量子化學(xué)計算需要進(jìn)行校正,沒有必要將同一全氮分子內(nèi)的氮原子類型進(jìn)行精細(xì)的分類。因此,全氮分子的范德華參數(shù)優(yōu)化均采用如下設(shè)置: 無定形結(jié)構(gòu)構(gòu)建時盒子密度設(shè)置為0.6 g·cm-3,對二聚體結(jié)合能進(jìn)行校正,每個全氮分子只定義一種原子類型。

3.4 參數(shù)特異性

表3列出了籠形全氮分子經(jīng)優(yōu)化后的范德華參數(shù),與Dreiding力場相比,優(yōu)化的范德華參數(shù)RMSE降低了約一個數(shù)量級,準(zhǔn)確程度大幅提高。ε集中在0.4～0.7 kJ·mol-1,σ集中在3.2～3.4 ?,說明對于不同的籠形結(jié)構(gòu),全氮分子的范德華參數(shù)確實存在特異性,只采用一種普適的氮原子參數(shù)不能準(zhǔn)確地描述所有的籠形全氮結(jié)構(gòu)。另外,隨著氮原子數(shù)的增加,ε與σ并沒有呈現(xiàn)明顯的規(guī)律。其中,N4(Td)的ε最大,σ也最大; N8(Oh)的ε最小,σ僅大于N12(D6h); N12(D6h)的σ最小,ε僅小于N4(Td)。究其原因,是籠形全氮分子特殊的空間結(jié)構(gòu)所致。

為了解釋籠形全氮分子范德華參數(shù)的特異性,采用MP2/cc-pvtz方法分別計算了籠形全氮分子棱-棱、底面-底面和側(cè)面-側(cè)面三種空間取向下的二聚體結(jié)合能與相對距離的關(guān)系圖,結(jié)果示于圖2。由圖中可知,不同空間取向下籠形全氮分子的平衡距離是較為接近的,差異主要體現(xiàn)在勢阱深度上。對于棱-棱取向,N12(D6h)勢阱深度最大,其次為N10(D5h)與N4(Td); 對于底面-底面取向,盡管N4(Td)底面上只有3個氮原子,卻顯示出了最強(qiáng)的分子間相互作用; 對于側(cè)面-側(cè)面取向,除了N4(Td)外,其余籠形結(jié)構(gòu)側(cè)面上均有4個氮原子,N12(D6h)結(jié)合能最大,N10(D5h)與N4(Td)次之。綜合所有取向可知,N4(Td)與N12(D6h)作用勢最大,這與其空間結(jié)構(gòu)特點是密切相關(guān)的,同時也可看出優(yōu)化的范德華參數(shù)是合理的。

表3籠形全氮分子范德華參數(shù)

Table3The van der Waals parameters of all-nitrogen molecules with cage structure

moleculesamplenumbervdWparametersε/kJ·mol-1σ/?RMSEDreidingthisworkN4(Td)45070．67993．37620．2110．034N6(D3h)32750．59963．36660．1810．099N8(Oh)30630．47153．21721．3320．372N10(D5h)31330．50563．24431．1190．386N12(D6h)33100．61893．21051．2850．395

a. edge-edge orientation

b. undersurface-undersurface orientation

c. side face-side face orientation

圖2不同空間取向的籠形全氮二聚體結(jié)合能與相對距離關(guān)系圖

Fig.2The binding energies of all-nitrogen dimer with cage structure as a function of distance for different orientations

3.5 晶體密度

基于建立的分子力場參數(shù),預(yù)測了20種全氮分子的晶體密度,為了便于和英國QinetiQ公司、瑞典國防研究院FOI的計算結(jié)果相比較,只列出了N4(Td)、N6(D3h)、N8(Oh)、N10(D5h)及N12(D6h)共計5種籠形全氮分子的晶體密度,結(jié)果依次為1.81,2.08,2.47,2.46,2.57 g·cm-3。除N4(Td)外,其余全氮分子的晶體密度均在2.0 g·cm-3以上,其中N12(D6h)晶體密度最高。由圖3可見,本研究結(jié)果介于文獻(xiàn)計算結(jié)果的范圍之內(nèi)。隨著氮原子數(shù)的增加,本研究預(yù)測的晶體密度并沒有呈簡單的增長趨勢,而是在N8(Oh)出現(xiàn)了一個突變點,這表明特異性力場參數(shù)的計算結(jié)果更加合理。

圖3晶體密度預(yù)測值與籠形全氮原子數(shù)關(guān)系圖

Fig.3Predicted crystal densities of all-nitrogen molecules with cage type as a function of nitrogen atom number

4 結(jié) 論

(1) 采用量子化學(xué)計算方法,建立了適用于全氮結(jié)構(gòu)的“特異性”力場參數(shù),即針對特定的全氮分子,建立了一套獨立的力場參數(shù)。由5種籠形全氮結(jié)構(gòu)棱-棱、底面-底面、側(cè)面-側(cè)面取向的二聚體結(jié)合能可知,不同結(jié)構(gòu)的全氮分子間相互作用確實存在差異,如只采用一套通用的力場參數(shù)勢必帶來較大誤差。因此,特異性力場參數(shù)在原理上更加精確合理。

(2) 范德華力是全氮分子在晶體堆積中的主導(dǎo)作用力,本研究詳細(xì)討論了范德華參數(shù)的優(yōu)化策略,包括二聚體抽樣個數(shù),二聚體結(jié)合能是否校正,參數(shù)種類等影響因素,與Dreiding力場中范德華參數(shù)相比,通過特異性參數(shù)計算的二聚體結(jié)合能與量子化學(xué)計算結(jié)果的均方根誤差降低了約一個數(shù)量級。

(3) 借助特異性力場參數(shù),預(yù)測了20種全氮分子的晶體密度。其中,N4(Td)、N6(D3h)、N8(Oh)、N10(D5h)及N12(D6h)五種籠形全氮分子的晶體密度分別為1.81，2.08，2.47，2.46,2.57 g·cm-3,介于英國QinetiQ公與瑞典國防研究院FOI報道結(jié)果的范圍之間,體現(xiàn)了特異性力場參數(shù)的真實性?；诒狙芯款A(yù)測的晶體密度,有助于進(jìn)一步探索籠形全氮分子的爆轟性能。

參考文獻(xiàn):

[1] Eremets M I, Gavriliuk A G, Trojan I A, et al. Single-bonded cubic form of nitrogen[J].NatureMaterials, 2004, 3(8): 558-563.

[2] Samartzis P C, Wodtke A M. All-nitrogen chemistry: how far are we from N60?[J].InternationalReviewsinPhysicalChemistry, 2006, 25(4): 527-552.

[3] Hirshberg B, Gerber R B, Krylov A I. Calculations predict a stable molecular crystal of N8[J].NatureChemistry, 2014, 6(1): 52-56.

[4] 李玉川, 龐思平. 全氮型超高能含能材料研究進(jìn)展[J]. 火炸藥學(xué)報, 2012, 35(1): 1-8.

LI Yu-chuan, PANG Si-ping. Progress of all-nitrogen ultrahigh-energetic material[J].ChineseJournalofExplosive&Propellants, 2012, 35(1): 1-8.

[5] 張光全, 董海山. 氮簇合物——潛在的高能量密度材料候選物[J]. 含能材料, 2004, 12(增刊): 105-113.

ZHANG Guang-quan, DONG Hai-shan. Nitrogen clusters—potential candidates as high-energy density materials[J].ChineseJournalofEnergeticMaterials(HannengCailiao), 2004, 12(Suppl.): 105-113.

[6] Tarver C M. Density estimations for explosives and related compounds using the group additivity approach[J].JournalofChemicalandEngineeringData, 1979, 24(2): 136-145.

[7] Ammon H L, Mitchell S. A new atom/functional group volume additivity data base for the calculation of the crystal densities of C, H, N, O and F-containing compounds[J].Propellants,Explosives,Pyrotechnics, 1998, 23(5): 260-265.

[8] Piacenza G, Legsa? G, Blaive B, et al. Molecular volumes and densities of liquids and solids by molecular mechanics—estimation and analysis[J].JournalofPhysicalOrganicChemistry, 1996, 9(6): 427-432.

[9] Wong M W, Wiberg K B, Frisch M J. Ab initio calculation of molar volumes: comparison with experiment and use in solvation models[J].JournalofComputationalChemistry, 1995, 16(3): 385-394.

[10] Politzer P, Martinez J, Murray J S, et al. An electrostatic interaction correction for improved crystal density prediction[J].MolecularPhysics, 2009, 107(19): 2095-2101.

[11] Holden J R, Du Z Y, Ammon H L. Prediction of possible crystal structures for C—, H—, N—, O—, and F—containing organic compounds[J].JournalofComputationalChemistry, 1993, 14(4): 422-437.

[12] Dzyabchenko A V, Pivina T S, Arnautova E A. Prediction of structure and density for organic nitramines[J].JournalofMolecularStructure, 1996, 378(2): 67-82.

[13] Haskins P J, Fellows J, Cook M D, et al. Molecular level studies of polynitrogen explosives[C]∥12thInternational Detonation Symposium, California, 2002.

[14] ?stmark H. High energy density materials (HEDM): overview, theory and synthetic efforts at FOI[C]∥New Trends in Research of Energetic Materials, Czech Republic, 2006: 231-250.

[15] Lee C, Yang W, Parr R G. Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density[J].PhysicalReviewB, 1988, 37(2): 785-789.

[16] Becke A D. Density-functional thermochemistry. Ⅲ. The role of exact exchange[J].JournalofChemicalPhysics, 1993, 98(7): 5648-5652.

[17] Frisch M J, Trucks G W, Schlegel H B, et al. Gaussian 09[C]∥Gaussian, Inc., Wallingford CT, 2009.

[18] Humphrey W, Dalke A, Schulten K. VMD: visual molecular dynamics[J].JournalofMolecularGraphics, 1996, 14(1): 33-38.

[19] Kramer C, Gedeck P, Meuwly M. Multipole-based force fields from ab initio interaction energies and the need for jointly refitting all intermolecular parameters[J].JournalofChemicalTheoryandComputation, 2013, 9(3): 1499-1511.

[20] Material Studio 8.0[C]∥Acceryls Inc.: San Diego, 2014.

[21] Mayo S L, Olafson B D, Goddard W A. Dreiding: A generic force field for molecular simulations[J].JournalofPhysicalChemistry, 1990, 94(26): 8897-8909.

[22] Head-Gordon M, Pople J A, Frisch M J. MP2 energy evaluation by direct methods[J].ChemicalPhysicsLetters, 1988, 153(6): 503-506.

[23] Kendall R A, Dunning Jr. T H, Harrison R J. Electron affinities of the first-row atoms revisited. Systematic basis sets and wave functions[J].JournalofChemicalPhysics, 1992, 96(9): 6796-6806.

[24] Boys S F, Bernardi F. Calculation of small molecular interactions by differences of separate total energies-some procedures with reduced errors[J].MolecularPhysics, 1970, 19(4): 553-566.