TATB晶體的狀態(tài)方程與振動性質(zhì)的密度泛函理論研究

蔣文燦， 陳 華， 張偉斌

(1. 中國工程物理研究院化工材料研究所，四川 綿陽 621999； 2. 中國工程物理研究院研究生部，四川 綿陽 621999)

1 引 言

a. crystal structure b. molecule structure

圖1 TATB晶體與分子結(jié)構(gòu)

Fig.1 Crystal and molecule structure for TATB

2 計算方法

理論計算運用VASP軟件[26]進(jìn)行，初始晶體結(jié)構(gòu)XRD實驗數(shù)據(jù)[1]，基于廣義梯度近似(GGA)的投影綴加平面波(Projector augmented wave，PAW)贗勢，引入vdW-DF2修正范德華力。截斷能設(shè)置為600 ev，k點取為2×2×2，截斷能和k點選擇保證每個原子最大受力變化不超過0.01 eV/?，結(jié)構(gòu)優(yōu)化以及狀態(tài)方程計算電子弛豫的標(biāo)準(zhǔn)為1×10-6eV/atom，原子弛豫標(biāo)準(zhǔn)1×10-5eV/atom，振動頻率通過DFPT方法[27]計算，收斂標(biāo)準(zhǔn)進(jìn)一步提高，電子弛豫標(biāo)準(zhǔn)為1×10-8eV/atom，原子最大受力不超過0.001 eV/ ?。

3 結(jié)果與討論

3.1 TATB晶體狀態(tài)方程

TATB晶體狀態(tài)方程(p-V曲線)擬合考慮兩種半經(jīng)驗狀態(tài)方程，即Vinet方程[28]和三階方程[29]，狀態(tài)方程形式如下:

(1)

(2)

a. equation of state for TATB crystal

b. cell length of a axis and c axisas a function of pressure

c. the ratio of cell parameters of b axisto a axis as a function of pressure

d. cell angle for α、β and γ as a function of pressure

圖2 本文計算的TATB晶體狀態(tài)方程以及晶胞參數(shù)與Olinger和Cady[30]以及Stevens等[4]實驗數(shù)據(jù)的對比

Fig.2 Comparison of the equation of state and cell parameters for TATB determined by Olinger and Cady[30], Stevens et al[4]and calculated by this work

圖3 TATB分子二面角隨壓力變化曲線

Fig.3 Dihedral angle of TATB molecular as a function of pressure

methodpressurerange/GPaBirch?MurnaghanK0K′0VinetK0K′0DFT?D2[31]0-818．77．9--DFT?D2[22]0-2013．611．715．19．2DFT?D3(BJ)[22]0-2013．88．714．37．9MD[32]0-1011．117．2--SAPT(DFT)[33]0-1012．411．0--OlingerandCady[30]0-7．016．75．716．65．9Stevensetal．[4]0-13．217．18．117．57．6Stevensetal．[4]0-5．013．612．414．710．0thiswork0-8．517．78．117．97．8

由圖2d可見，加壓過程中TATB晶體β和γ晶胞角幾乎不變，由圖2b可見c軸晶胞參數(shù)壓縮性高于a軸，圖2c可見a軸和b軸晶胞參數(shù)壓縮率為常數(shù)，這與Olinger和Candy[30]以及其他理論計算結(jié)果[22,31]相符，表明Olinger和Candy[30]提出的兩個假定是合理的。Ojede和Cagin[6]利用DFT計算結(jié)果發(fā)現(xiàn)加壓過程在1.5 GPa附近TATB晶體內(nèi)c軸方向氫鍵作用數(shù)明顯增加，本研究統(tǒng)計c軸方向相鄰層氫鍵作用數(shù)發(fā)現(xiàn)(截斷半徑2.5 ?)，8.5 GPa時氫原子與相鄰層氧原子最近距離由3.040 ?變化為2.686 ?，表明氫鍵接觸數(shù)量并沒有增加。但隨著壓力進(jìn)一步增大，c軸方向氫鍵作用數(shù)可能會增多。由圖3可見0～8.5 GPa加壓過程中苯環(huán)二面角C(1)—C(2)—C(3)—C(4)由3°變?yōu)?.7°，表明壓力將會引起碳原子垂直苯環(huán)平面的扭轉(zhuǎn)(即折疊)。另外，加壓過程中氧二面角O(1)—N(2)—C(2)—C(1)由0.4°變?yōu)?.9°，而氫原子二面角H(6)—C(5)—N(5)—C(6)變化較小，由0.7°變化為1.5°，表明壓力引起氧原子平面外扭轉(zhuǎn)。在2 GPa附近二面角出現(xiàn)拐點，拐點出現(xiàn)原因有待進(jìn)一步探索。以上計算結(jié)果表明，TATB在加壓過程中出現(xiàn)了兩處微結(jié)構(gòu)變化，使TATB分子相互彎曲靠近。

3.2 TATB晶體振動性質(zhì)

3.2.1 零壓TATB晶體振動性質(zhì)

TATB晶體中含有48個原子，根據(jù)群論相關(guān)知識[34]，共有144種振動模，其中72種拉曼振動模和72種紅外振動模，由于紅外光譜[7]能觀察到的振動峰較少，而拉曼實驗可以給出TATB晶體大部分振動峰[35]，且低頻段普遍存在晶格振動與分子振動耦合現(xiàn)象，因此對振動頻率高于250 cm-1的拉曼振動模式進(jìn)行了分配，振動模式指認(rèn)結(jié)果見附表(附表可查閱本刊網(wǎng)站)。表2列出了與Liu等[15]分配結(jié)果出現(xiàn)差異的振動模式以及下文高壓振動性質(zhì)研究中重點關(guān)注的振動模。分析表2發(fā)現(xiàn)，Q39和Q42處Liu等[15]對TATB分子指認(rèn)結(jié)果分別為C—NO2伸縮振動與NO2擺動(表2中未列出)，但對TATB晶體分配結(jié)果為NO2擺動和C—NO2伸縮振動(見表2)，振動形式完全相反，而本文指認(rèn)結(jié)果與Liu等[15]對Q39和Q42處TATB分子振動指認(rèn)結(jié)果一致。Q50處振動模式除Liu等[15]觀察到的NH2擺動以外，還出現(xiàn)了NO2擺動，因此振動情況更加復(fù)雜。Q65處振動模式與Liu等[15]結(jié)果完全不同，指認(rèn)結(jié)果為苯環(huán)扭動。Q70處與Liu等[15]振動分配出現(xiàn)區(qū)別，為苯環(huán)面外變形振動，而隨著振動頻率進(jìn)一步增加，Q72和Q74處苯環(huán)面外變形振動變?nèi)?，而NH2面外彎曲振動加強(qiáng)并逐漸變?yōu)橹饕駝有问剑擞嬎憬Y(jié)果與Liu等[15]不同，但與Sui等[10]紅外實驗對730 cm-1處振動分配結(jié)果相符(Q72)。Sui等[10]實驗觀察到783 cm-1附近振動主要為NO2剪切振動，Liu等[15]結(jié)果為NO2扭曲振動，而本文結(jié)果為NH2平面外彎曲振動和NO2剪切振動(Q80，Q81和Q84)，與Sui等[10]實驗結(jié)果符合更好。1120 cm-1處(Q100)Sui等[10]人紅外實驗結(jié)果表明為C—NO2伸縮振動，本文指認(rèn)結(jié)果為C—NO2伸縮振動，與Sui等[10]人實驗結(jié)果相符，而Liu等[15]指認(rèn)結(jié)果為NH2和NO2振動。另外，TATB晶體內(nèi)振動耦合情況明顯，如Q88，Q89，Q91，Q94，Q95指認(rèn)結(jié)果表明，Q88和Q89在NH2擺動和NO2剪切振動基礎(chǔ)上，可能還含有C—NO2伸縮，而隨著振動頻率增加，在Q91、Q94和Q95中出現(xiàn)了明顯的C—NO2伸縮振動，而苯環(huán)伸縮振動也逐漸變得明顯，逐漸變得復(fù)雜的振動耦合使振動模分配結(jié)果與Liu等[15]以及Sui等[10]實驗對振動模式分配結(jié)果出現(xiàn)差異。如1030 cm-1處(Q94)，本文分配結(jié)果為NH2擺動、NO2剪切振動以及C—NO2伸縮振動，可能還含有苯環(huán)伸縮振動，Liu等[15]結(jié)果為NH2擺動和C—NO2伸縮振動，Sui等[10]和Town等[12]均為苯環(huán)伸縮振動。更加復(fù)雜的振動模式出現(xiàn)在1086～1300 cm-1之間，主要表現(xiàn)為苯環(huán)振動與NH2和NO2振動耦合，如Q105處振動同時包含苯環(huán)伸縮振動、NH2剪切振動、NH2擺動以及C—NO2伸縮振動， Q107處振動模式指認(rèn)結(jié)果為苯環(huán)伸縮，C—NO2伸縮振動，C—NH2伸縮振動與擺動，可能還含有NO2的反對稱伸縮振動，而這種復(fù)雜振動也可能是多種基頻模發(fā)生費米共振引起[10]，隨著外界壓力增加，這種復(fù)雜的振動耦合很可能引起組分變化并導(dǎo)致振動峰劈裂與合并，并在加壓和減壓過程可能表現(xiàn)為不可逆變化。另外，3000～4000 cm-1振動頻率與實驗數(shù)據(jù)[35]相差較大(見附表)，最大偏差達(dá)到157 cm-1，而Liu等[15]DFT-GGA計算結(jié)果在3000～4000 cm-1之間與實驗相比也出現(xiàn)了130 cm-1左右的偏差，這種偏差可能與DFT方法有關(guān)。

表2 TATB晶體部分拉曼活性模振動頻率及其分配

Table 2 Part of Raman active vibration frequencies and their assignment in TATB crystal

modevibrationfrequency/cm-1Liu，etal[15]thisworkQ30291ringtwistringtwistQ36362ringdefromationringdefromationQ39381NO2rockC—NO2stretchQ42384C—NO2stretchNO2rockQ50522NH2rockNH2rock+NO2rockQ65700C—NO2stetchringtwistQ70711C—NO2torsionringdistortionQ72729C—NO2torsionNH2wagQ74753ringdistortion+NH2wagNH2wagQ76758C—NH2torsionringdistortion+NH2wagQ78791NH2wagNH2wagQ80794NH2twistNH2twist+NO2scissorQ81797C—NH2stretchNH2twist+NO2scissorQ84805NH2twistNH2twist+NO2scissorQ86817NH2twistNH2twistQ88849NH2rock+NO2scissorNH2rock+NO2scissor+C—NO2stretchQ89852NH2rock+NO2scissorNH2rock+NO2scissor+C—NO2stretchQ911000NH2rock+C—NO2stretchNH2rock+C—NO2stretch+NO2scissorQ941003NH2rock+C—NO2stretchNH2rock+C—NO2stretch+NO2scissorQ951087C—NO2stretchNH2rock+C—NO2stretch+NO2scissorQ1001120ringstretch+C—NH2stretchringstretch+C—NH2stretch+NH2scissorQ1071164ringstretchringstretch+C—NO2stretch+NH2rock+C—NH2stretch

3.2.2 加壓TATB晶體振動性質(zhì)

由于高壓拉曼實驗僅給出100～900 cm-1振動頻率隨壓力變化結(jié)果，又因為1000～2000 cm-1振動形式比較復(fù)雜，本文分析了290～900 cm-1加壓過程振動耦合情況，為了便于比較，使拉曼實驗[9]觀察到的振動頻率數(shù)目與理論計算結(jié)果數(shù)量一致，結(jié)果如圖4。對比發(fā)現(xiàn)，250～900 cm-1頻段理論計算振動頻率與實驗數(shù)據(jù)符合較好。為了探究壓力對TATB晶體振動頻率的影響，分析了Q30、Q36、Q50、Q65和Q81處振動模式同時對250～900 cm-1頻段其他振動頻率處振動模式進(jìn)行了指認(rèn)，對0～8.5 GPa加壓過程振動模式分析發(fā)現(xiàn)，隨著壓力增加，除Q39，Q42，Q86，Q76和Q78以外，其余振動頻率處振動形式基本保持不變，加壓過程在290～900 cm-1之間并未發(fā)現(xiàn)振動頻率突變。對Q39和Q42振動模式分析發(fā)現(xiàn)，加壓過程中在1.24 GPa時Q39由C—NO2伸縮振動變?yōu)镹O2擺動，在8.5 GPa時，Q39變?yōu)楸江h(huán)變形振動，Q42處振動模式在1.93 GPa由硝基擺動變?yōu)镃—NO2伸縮振動，之后隨著壓力增加，振動形式不變。需要注意的是，Q39和Q42處均出現(xiàn)了苯環(huán)振動，但這種振動可能由NH2伸縮振動或者NO2擺動引起，也可能在Q42處和Q39處苯環(huán)與其余振動耦合在一起。Q76處在零壓時為NH2面外擺動和面外彎曲振動與苯環(huán)面外扭曲振動，在1.42 GPa時苯環(huán)振動減弱，而在6.38 GPa時變?yōu)镹H2面外擺動和NH2面外彎曲振動，含有NO2剪切振動，但并不明顯，在8.5 GPa時NO2剪切振動變得明顯。Q78處(如圖5)在零壓時為NH2面外擺動，但在0.5 GPa時還出現(xiàn)了NH2面外彎曲振動和NO2剪切振動，但并不明顯，在1.24 GPa時NH2面外彎曲振動和NO2剪切振動變得明顯。Q86處振動模式變化比較明顯，在零壓時為NO2面外彎曲振動，在4.67 GPa時為NH2面外彎曲振動和NO2剪切振動耦合。分析以上振動模式變化情況發(fā)現(xiàn)，隨著壓力增加，NH2面外彎曲振動或擺動與NO2剪切振動耦合，而對TATB分子結(jié)構(gòu)分析結(jié)果表明在加壓過程中相鄰層分子相互彎曲靠近，表明TATB這種微結(jié)構(gòu)變化導(dǎo)致了振動耦合，分子間氫鍵作用增強(qiáng)。

圖4 加壓過程中TATB晶體振動頻率的計算結(jié)果(紅色圓點)與實驗結(jié)果[9](藍(lán)色線)的對比

Fig.4 Comparison of the calculated values (the red dot) and the experimental ones[9](the blue line) for the vibrational frequencies under pressure process of TATB crystal

a. 0 GPa b. 8.5 GPa

圖5 TATB晶體0 GPa和8.5 GPa時791 cm-1處振動模指認(rèn)結(jié)果(圖中紅色代表氧原子，藍(lán)色代表氮原子，白色代表氫原子，黑色代表碳原子)

Fig.5 The identification results of the vibration mode of TATB crystal at 0 GPa and 8.5 GPa at 791 cm-1. The red ball in picture represents oxygen atom. The blue ball represents nitrogen atom. The white ball represents hydrogen atom. The black ball represents carbon atom

4 結(jié) 論

利用vdW-DF2研究了TATB晶體狀態(tài)方程以及振動性質(zhì)，理論計算結(jié)果與實驗符合較好，表明vdW-DF2能夠較好地表征TATB晶體結(jié)構(gòu)與性質(zhì)。對狀態(tài)方程研究發(fā)現(xiàn)，加壓過程中TATB出現(xiàn)了兩處微結(jié)構(gòu)變化，即苯環(huán)折疊與硝基扭轉(zhuǎn)，隨著壓力增加兩者的變化程度增大，并在2GPa附近出現(xiàn)拐點。對TATB晶體部分分子內(nèi)振動模式進(jìn)行了重新分配，研究結(jié)果表明在1100～1500 cm-1波數(shù)之間TATB晶體振動尤其復(fù)雜，氨基與硝基和苯環(huán)振動耦合。結(jié)合理論擬合的TATB晶體狀態(tài)方程，研究了壓力為0～8.5 GPa，波數(shù)為290～900 cm-1，TATB晶體振動耦合情況以及分子間相互作用過程，研究發(fā)現(xiàn)，隨著壓力增加，TATB相鄰層分子相互彎曲靠近，引起氨基平面外彎曲振動或擺動與硝基剪切振動耦合，表明分子間氫鍵作用增強(qiáng)。

參考文獻(xiàn):

[1] Cady H H, Larson A C. The crystal structure of 1, 3, 5-triamino-2, 4, 6-trinitrobenzene[J].ActaCrystallographica, 1965, 18(3): 485-496.

[2] Phillips DS, Schwarz RB, Skidmore CB, et al. Some observations on the structure of TATB[C]∥Shock compression of condensed matter-1999. AIP Publishing, 2000: 707-710.

[3] Kolb J R, Rizzo H F. Growth of 1,3,5-triamino-2,4,6-trinitrobenzene (TATB) I. Anisotropic thermal expansion[J].Propellants,Explosives,Pyrotechnics, 1979, 4(1): 10-16.

[4] Stevens L L, Velisavljevic N, Hooks D E, et al. Hydrostatic compression curve for triamino-trinitrobenzene determined to 13.0 GPa with powder X-ray diffraction[J].Propellants,Explosives,Pyrotechnics, 2008, 33(4): 286-295.

[5] Manaa M R, Fried L E. Nearly equivalent inter- and intramolecular hydrogen bonding in 1,3,5-triamino-2,4,6-trinitrobenzene at high pressure[J].TheJournalofPhysicalChemistryC, 2011, 116(3): 2116-2122.

[7] Pravica M, Yulga B, Liu Z, et al. Infrared study of 1, 3, 5-triamino-2, 4, 6-trinitrobenzene under high pressure[J].PhysicalReviewB, 2007, 76(6): 64102.

[8] Pravica M, Yulga B, Tkachev S, et al. High-pressure far-and mid-infrared study of 1, 3, 5-triamino-2, 4, 6-trinitrobenzene[J].TheJournalofPhysicalChemistryA, 2009, 113(32): 9133-9137.

[9] Davidson A J, Dias R P, Dattelbaum D M, et al. “stubborn” Triaminotrinitrobenzene: unusually high chemical stability of a molecular solid to 150 GPa[J].TheJournalofChemicalPhysics, 2011, 135(17): 174507.

[10] Sui H, Zhong F, Cheng K, et al. IR vibrational assignments for 1, 3, 5-triamine-2, 4, 6-trinitrobenzene (TATB) based on the temperature-dependent frequency shifts[J].SpectrochimicaActaPartA:MolecularandBiomolecularSpectroscopy, 2013, 114(10): 137-143.

[11] 劉紅，趙紀(jì)軍，龔自正，等. 原子與分子物理學(xué)報，2007, 24 (2), 291-297.

LIU Hong, ZHAO Ji-jun, GONG Zi-zheng, et al. Crystal structure of the TATB crystal under high pressure[J].JournalofAtomicandMolecualrPhysics, 2007, 24(2): 291-297.

[12] Towns TG. Vibrational spectrum of 1,3,5-triamino-2,4,6-trinitrobenzene[J].SpectrochimicaActaPartaMolecularSpectroscopy, 1983, 39(9): 801-804.

[13] Sorescu DC, Rice BM. Theoretical predictions of energetic molecular crystals at ambient and hydrostatic compression conditions using dispersion corrections to conventional density functionals (DFT-D)[J].TheJournalofPhysicalChemistryC, 2010, 114(14): 6734-6748.

[14] 劉紅. 含能材料的高壓行為和光譜學(xué)特性的原子模擬[D] . 四川: 西南交通大學(xué)，2006.

[15] Liu H, Zhao J, Ji G, et al. Vibrational properties of molecule and crystal of TATB: a comparative density functional study[J].PhysicsLettersA, 2006, 358(1): 63-69.

[17] Grimme S. Semiempirical GGA-type density functional constructed with a long-range dispersion correction[J].JournalofComputationalChemistry, 2006, 27(15): 1787-1799.

[18] Lee K, Murray é D, Kong L, et al. Higher-accuracy van der waals density functional[J].PhysicalReviewB, 2010, 82(8): 81101.

[19] Grimme S, Antony J, Ehrlich S, et al. A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu[J].TheJournalofChemicalPhysics, 2010, 132(15): 154104.

[20] Dion M, Rydberg H, Schr?der E, et al. Van der waals density functional for general geometries[J].PhysicalReviewLetters, 2004, 92(24): 246401.

[21] Grimme S, Ehrlich S, Goerigk L. Effect of the damping function in dispersion corrected density functional theory[J].JournalofComputationalChemistry, 2011, 32(7): 1456-1465.

[22] Fedorov IA, Zhuravlev YN. Hydrostatic pressure effects on structural and electronic properties of TATB from first principles calculations[J].ChemicalPhysics, 2014, 436(7): 1-7.

[23] Shimojo F, Wu Z, Nakano A, et al. Density functional study of 1, 3, 5-trinitro-1, 3, 5-triazine molecular crystal with van der waals interactions[J].TheJournalofChemicalPhysics, 2010, 132(9): 94106.

[24] Berland K, Hyldgaard P. Analysis of van der waals density functional components: binding and corrugation of benzene and C 60 on boron nitride and graphene[J].PhysicalReviewB, 2013, 87(20): 205421.

[25] Dobson JF, Gould T. Calculation of dispersion energies[J].JPhys:CondensMatter, 2012, 24(7): 73201.

[26] Kresse G, Furthmüller J. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set[J].PhysicalReviewB, 1996, 54(16): 11169.

[27] Baroni S, De gironcoli S, Dal corso A, et al. Phonons and related crystal properties from density-functional perturbation theory[J].ReviewsofModernPhysics, 2001, 73(2): 515.

[28] Vinet P, Ferrante J, Smith J, et al. A universal equation of state for solids[J].JournalofPhysicsC:SolidStatePhysics, 1986, 19(20): L467.

[29] Birch F. Finite elastic strain of cubic crystals[J].PhysicalReview, 1947, 71(11): 809.

[30] Olinger BW, Cady HH. Hydrostatic compression of explosives and detonation products to 10 GPa (100 Kbars) and their calculated shock compression: results for PETN, TATB, CO2, and H2O[R]. Los Alamos Scientific Lab., N. Mex.(USA), 1976.

[31] Budzevich M, Landerville A, Conroy M, et al. Hydrostatic and uniaxial compression studies of 1, 3, 5-triamino-2, 4, 6-trinitrobenzene using density functional theory with van der waals correction[J].JournalofAppliedPhysics, 2010, 107(11): 113524.

[32] 金釗, 劉建, 王麗莉, 等. 物理化學(xué)學(xué)報，2014, 30 (4): 654-661.

JIN Zhao, LIU Jian, WANG Li-li, et al. Development and validation of an all-atom force field for the energetic materials TATB, RDX and HMX[J].ActaPhysico-chimicaSinica, 2014, 30(4): 654-661.

[33] Taylor DE. Intermolecular forces and molecular dynamics simulation of 1,3,5-triamino-2,4,6-trinitrobenzene (TATB) using symmetry adapted perturbation theory[J].TheJournalofPhysicalChemistryA, 2013, 117(16): 3507-3520.

[34] 吳國禎. 分子振動光譜學(xué)[M]. 北京: 清華大學(xué)出版社，2001: 98-110.

[35] Mcgrane S, Shreve A. Temperature-dependent raman spectra of triaminotrinitrobenzene: anharmonic mode couplings in an energetic material[J].TheJournalofChemicalPhysics, 2003, 119(12): 5834-5841.