含能材料的量子化學計算與分子動力學模擬綜述

居學海，葉財超，徐司雨

（1.南京理工大學化工學院，江蘇 南京210094；2.西安近代化學研究所燃燒與爆炸技術重點實驗室，陜西 西安710065）

引 言

近30年來，量子化學計算方法和分子模擬技術、以及計算機技術的飛速發(fā)展，對材料科學的發(fā)展產(chǎn)生了深刻影響。材料科學與物理、化學、數(shù)學和工程力學等學科相互交叉與滲透，產(chǎn)生了計算材料學這一新興學科。通過對材料組分和材料結構的模型化計算實現(xiàn)對材料設計、制備、加工、結構和性能等參數(shù)或過程的定量描述，建立材料結構與性能之間的關系，并最終按指定目標設計新材料。研究對象的空間和時間尺度不同，計算材料學的研究范圍、深度和方法也隨之不同。按材料的空間尺度大體可劃分為4個層次：電子／原子層次、微觀、介觀和連續(xù)體（宏觀）層次；或簡單地分為微觀、介觀和宏觀3個層次。計算方法依次為量子力學方法（半經(jīng)驗分子軌道法、第一性原理）、第一性原理分子動力學模擬（QMD）、經(jīng)典分子動力學（MD）和蒙特卡羅模擬、耗散粒子動力學模擬、相場動力學和原胞自動化方法、有限元和有限差分法等。一般來說，通過電子及原子層次計算實現(xiàn)對材料本質的認識，是材料設計的科學基礎。其研究方法主要為量子力學（在化學和材料學領域，常稱量子化學）和分子動力學模擬。量子化學（QC）的研究對象為平衡態(tài)、單分子或幾個分子組成的體系、周期性體系（晶胞或超晶胞）；不適用于動力學過程和有溫度壓力變化的體系。分子動力學模擬描述原子核的運動過程，求得系統(tǒng)的結構和性質，模擬結果既包含系統(tǒng)的靜態(tài)性質，又包含動態(tài)特性；經(jīng)典（經(jīng)驗勢）分子動力學可研究上百萬甚至千萬個原子的體系，但不適用于有電子轉移、原子變價的過程。后者須采用第一性原理分子動力學。本文主要介紹量子化學和經(jīng)典分子動力學模擬在含能材料研究中的應用及進展。

1 熱力學性質及爆轟性能

含能材料的爆轟性能與爆熱密切相關，而爆轟過程的反應熱可由爆轟反應各物質的生成熱（HOFs）直接求得。通過選用適當?shù)姆椒?，準確計算含能材料分子的生成熱是量子化學的一大優(yōu)勢。最早使用半經(jīng)驗分子軌道方法可迅速并直接計算出生成熱。但這些半經(jīng)驗方法主要依據(jù)一些代表性小分子和烴類的熱力學和光譜數(shù)據(jù)進行參數(shù)化，對含多種取代基或特殊結構的含能材料生成熱的計算往往帶來較大的誤差。一些改進的半經(jīng)驗方法，如PDDG（Pairwise Distance Directed Gaussian modification）改進法，可減少誤差。比較了622種含CHNO化合物生成熱計算值與實驗值的絕對誤差，PM3方法為18.4kJ／mol，PDDG／PM3方法為13.4kJ／mol［1］。

隨著計算技術的發(fā)展，第一性原理計算（從頭計算和密度泛函）逐漸取代半經(jīng)驗分子軌道法在含能材料中普遍運用。但第一性原理只能求得分子的總能量，無法直接計算生成熱。這就要求設計等鍵反應，利用參考物的實驗生成熱，借助Hess定律，求得目標分子的生成熱。在設計等鍵反應時，根據(jù)鍵分離規(guī)則把分子分解成一系列與所求物質具有相同化學鍵類型的小分子（已知生成熱）。為減少誤差，設計等鍵反應時應盡量使母體骨架或原有化學鍵保持不變。

第一性原理計算方法很多。含能材料閉殼層限制性Hartree-Fock（HF）計算可得到準確的分子幾何構型，但求得的能量與實驗誤差通常在200kJ／mol以內，個別誤差高達700kJ／mol［2］。因HF 方法忽略電子相關效應，導致能量值出現(xiàn)系統(tǒng)正誤差。為此，可用微擾法（MP2和MP4等）、多組態(tài)法（如CISD、CCSD 和CASSCF）及密度泛函理論（DFT）法（如B3LYP、B3PW91）校正電子相關效應。微擾法和DFT 法計算的能量誤差通常在40kJ／mol以內，而多組態(tài)法計算的誤差一般小于5～8kJ／mol［2］。多組態(tài)法即使對中等體系的計算量也非常大，因此通常用DFT 法處理電子相關問題。雖然DFT 法的能量絕對誤差40kJ／mol是一個不小的值，但如前所述，生成熱的計算通過設計等鍵反應來實現(xiàn)，反應物與產(chǎn)物的能量誤差大部分相互抵消。DFT 泛函的選用，通常要以同類型化合物的準確實驗生成熱為基準，比較并檢驗各種泛函對該類化學物的計算結果，以確定合適方案。

量子化學計算出的生成熱是氣相生成熱，結合遺傳算法（GA）、靜電勢法和神經(jīng)網(wǎng)絡法（NN）等求得固體的升華熱后，方可求得固相生成熱。靜電勢法的計算誤差小，在含能材料升華熱計算中應用較廣［3］。對有機分子組成的固體，其升華熱取決于分子間相互作用能，相互作用能越大則升華熱越大。因此，升華熱與分子表面的靜電勢有關。值得一提的是，新近基于1 384個化合物（含172個化學基團）升華熱所建立的3層前饋神經(jīng)網(wǎng)絡計算模型［4］，其升華熱的計算值與實驗結果的相關系數(shù)平方、平均誤差和均方根誤差分別為0.985 4，3.54%和4.21kJ／mol。遠小于靜電勢法的誤差11.7kJ／mol，具有一定的適用性和精確性。對于含能材料，該方法能否準確預測化合物的升華熱還有待檢驗。

用Kamlet經(jīng)驗公式［5］預測爆速和爆壓時，不僅需要生成熱值，還要有密度值。基于化合物優(yōu)化構型，用Monte-Carlo 法求得分子周邊0.001e／bohr3等電子密度面所包圍的體積，得到摩爾體積（V）。因該法求得的體積值在較大范圍內波動，通常要取100次以上平均值。由相對分子質量與平均摩爾體積之比求得理論密度（ρ）。計算表明，在B3LYP／6-31G＊＊水平下計算的理論密度與實測值吻合較好［6］。表1列出一些含能材料分子的氣相生成熱、爆速和爆壓。結果表明，計算值與實驗值相吻合。用第一性原理計算結合等鍵反應不僅求得生成熱，還能求得基團相互作用能和環(huán)張力［7-8］，從分子水平闡明結構與性能的關系。

表1 一些含能材料的氣相生成熱（ΔfHgas）、爆壓和爆速Table 1 Heats of formation，detonation pressure and velocity for some energetic materials

2 分解反應機理

含能材料的爆燃過程是以凝聚相和氣相化學反應（含氣固界面反應）為基礎的多組分多階段過程。量子化學計算是通過研究熱分解過程各基元反應歷程，建立反應物、過渡態(tài)（包括中間體）和產(chǎn)物之間的聯(lián)系。近年來，第一性原理方法廣泛應用于研究爆燃條件下的熱分解基元反應。例如，量子化學計算揭示氨基四唑的分解途徑主要有兩種方式［10］：兩鍵斷裂開環(huán)分解成RN3和NH2CN 及一鍵斷裂開環(huán)生成疊氮化物，該疊氮化物再失去一分子N2；前一種反應能壘較低。5-硝基-1-氫-四唑衍生物的熱分解包括直接開環(huán)途徑和質子轉移途徑［11］；開環(huán)可以是由N1-N2鍵斷裂開環(huán)，也可以是C5-N4鍵斷裂開環(huán)；開環(huán)之后，兩個途徑都由另一個N-N 鍵斷裂產(chǎn)生一個N2。密度泛函計算表明［12］，F(xiàn)OX-7不是通過C-NO2鍵直接斷裂分解，而是通過硝基-亞硝基重排機理分解。硝基-亞硝基重排反應的活化能計算值為247kJ／mol，接近FOX-7分解活化能的實測值242kJ／mol。這些研究表明，含能材料化學鍵斷裂生成小分子碎片之前，有時會發(fā)生氫遷移、異構化或重排反應。

用量子化學研究含能材料分解機理時，通常只考慮單分子過程，即將含能材料當作氣相體系。顯然，凝聚態(tài)中的鄰近分子間相互作用在很大程度上影響甚至改變含能材料的分解歷程。DFT 計算氣相TNAZ（1，3，3-三硝基氮雜環(huán)丁烷）分解時，其NNO2和C-NO2鍵斷裂能分別為159kJ／mol和171kJ／mol，接近氣相TNAZ分解活化能實驗值164.7kJ／mol；其HONO消除反應活化能分別為183.9和188kJ／mol，與固相分解活化能194.8kJ／mol相近；由此推測氣相分解主要是N-NO2和C-NO2鍵的斷裂，而固相分解主要是HONO 的消除［13］。隨后的計算進一步證實，TNAZ初始分解有NO2斷裂、HONO消去和開環(huán)3個過程。NO2斷裂的速率最大，但在固體中后兩個過程也有重要作用，即固體中鄰近分子間相互作用影響反應歷程［14］。用DFT 計算NTO 二聚體分解過程得知，二聚體經(jīng)歷活化能為367kJ／mol的基元反應之后，在總反應的第2步即生成CO2［15］；而在單分子分解反應中CO是終產(chǎn)物之一［16］。類似地，用DFT-B3LYP方法研究5-氨基四唑（5-ATZ）的分解過程時發(fā)現(xiàn)，雖然N2消去反應是單分子分解的主要反應途徑，其氫鍵二聚體能顯著地降低HN3消去反應的活化能，使HN3的消去反應優(yōu)先進行，這與實驗結果相一致［17］。

含能材料分解時產(chǎn)生高溫，導致激發(fā)態(tài)反應的發(fā)生，有時甚至占主導地位。探討激發(fā)態(tài)反應過程可揭示一些新的實驗現(xiàn)象。例如，NO是RDX 激發(fā)態(tài)（經(jīng)226～258nm 光子激發(fā)）分解時唯一觀測到的小分子產(chǎn)物［18］。S1和S2態(tài)垂直激發(fā)能的計算值分別為5.59和6.1eV，表明RDX分子吸收226～258nm（5.5～4.8eV）光子后被激發(fā)至第一激發(fā)態(tài)。計算結果表明［19］，其激發(fā)態(tài)分解須經(jīng)歷S1／S0圓錐交叉點（CI）的nitro-nitrite異構化非絕熱過程。第2 垂直激發(fā)態(tài)S2，F(xiàn)C分解時，必然經(jīng)歷S2／S1圓錐交叉點（S2／S1）CI。即從S2態(tài)分解時，其最小能壘通道為S2（FC）→（S2／S1）CI→（S1／S0）CI→S0（nitrite）→NO 消去。此外，還發(fā)現(xiàn)2個協(xié)同的nitro-nitrite異構化過渡態(tài)，表明分解同時產(chǎn)生3個NO分子，與硝胺類分子激發(fā)態(tài)分解時NO信號強度正比于N-NO2數(shù)量的實驗結果相符。圖1為RDX激發(fā)態(tài)分解所產(chǎn)生的NO發(fā)生［A2Σ+ （ν′=0）←X2Π（ν″=0，1，2，3）］躍遷時的REMPI光譜，其轉動和振動溫度分別為20K 和1 800K［20］。計算結果與REMPI 光譜相吻合。高振動溫度（1 800K）產(chǎn)物NO 的形成，表明激發(fā)態(tài)RDX 經(jīng)勢能面交叉點快速轉變成基態(tài)，并沿基態(tài)勢能面分解。HMX和CL-20的激發(fā)態(tài)分解與RDX 相似，也產(chǎn)生低轉動溫度及高振動溫度的NO。與此相反，實驗發(fā)現(xiàn)二甲硝胺經(jīng)226nm 光子激發(fā)并分解后，產(chǎn)生高轉動溫度（120K）NO，CASSCF的計算結果預測了上述現(xiàn)象［20］。

圖1 RDX 激發(fā)態(tài)分解產(chǎn)物NO 的高分辨REMPI譜（轉動和振動溫度分別為20和1 800K）Fig.1 REMPI spectra of the NO product from excited electronic state decomposition of RDX（Spectra rotational and vibrational temperature are 20and 1800K，respectively）

需要指出的是，由于目前無法對含能材料反應勢能面進行掃描，即無法確定勢能面上所有的鞍點和能量極小點，導致不同文獻對相同分子的分解機理報道不完全一致。當然，在少數(shù)情況下計算方法也導致結果的差別。故用量子化學方法探討反應機理時，應注重不同計算方法之間的相互比較和印證。

含能材料的分解過程涉及溫度壓力變化和多分子多組分碰撞。量子化學無法完整揭示這些過程，需要借助分子動力學模擬進行補充和完善。后者不僅揭示溫度和壓力的影響，還給出中間產(chǎn)物、后續(xù)反應過程和各物種隨時間的分布。對RDX 在高沖壓下的分子動力學模擬表明，低壓縮比時RDX較少分解，當壓縮比大于40%時RDX 快速分解，并預測了中間產(chǎn)物（NO2，NO，HONO，OH）和最終產(chǎn)物（H2O，N2，CO，CO2）的分布規(guī)律［21］?；诹孔恿W多尺度分子反應動力學模擬研究TATB受到9 km／s沖擊作用0.43ns及10km／s沖擊作用0.2ns時，發(fā)現(xiàn)TATB形成富氮雜環(huán)分子簇，從而阻礙了TATB進一步分解［22］。

QC與MD 相結合，可揭示起爆或燃燒過程中的反應歷程、決速步驟及活化能、催化機理、主次反應、產(chǎn)物分布、溫度和壓力等外界因素對反應過程的影響等。今后反應機理理論研究將更注重多分子和多組分體系、激發(fā)態(tài)過程以及界面反應。

3 感 度

20世紀70年代末，人們依據(jù)構效關系模型及其主要建模參數(shù)，相繼提出以下撞擊感度判據(jù)：氧平衡、激發(fā)態(tài)電荷梯度（ΔC＊／l）、最小鍵級、靜電勢、活化能、鍵離解能和化學位移等。除早期的氧平衡判據(jù)外，其他感度判據(jù)均為量子化學參數(shù)，并且主要針對撞擊感度而言。

早期炸藥感度研究中，Kamlet等歸納了大量有機炸藥實驗數(shù)據(jù)，提出感度-氧平衡理論。經(jīng)實驗數(shù)據(jù)擬合，建立了根據(jù)炸藥氧平衡計算撞擊感度的公式［23］。顯然，這一感度判據(jù)沒有直接涉及含能材料的化學屬性，特別是沒有涉及其電子結構，無法揭示感度的本質。

隨著理論研究的不斷發(fā)展，發(fā)現(xiàn)僅考慮炸藥分子內基團的特性和數(shù)目是不夠的，必須深入到電子結構。Depluech等用半經(jīng)驗量子化學方法計算RNO2類典型炸藥，發(fā)現(xiàn)基態(tài)電子結構和吸收能分布決定炸藥的起爆性能［24］，而這兩個因素又影響到激發(fā)態(tài)的電子結構，并建議用ΔC＊／l值（ΔC＊是激發(fā)態(tài)下R-NO2鍵兩端原子上凈電荷之差，l為RNO2鍵的鍵長；兩者的比值為激發(fā)態(tài)R-NO2鍵上的電荷梯度）作為硝基類炸藥的撞擊感度判據(jù)。這是從電子結構層次闡明炸藥爆炸性能的開拓性工作。但對結構差別較大的炸藥，其ΔC＊／l值與感度的相關性不顯著；另外并非任何外界作用都會使炸藥的分子從基態(tài)躍遷到激發(fā)態(tài)。將電子激發(fā)過程的難易程度用于判斷含能材料晶體的感度，發(fā)現(xiàn)撞擊感度與前沿帶隙相關，帶隙越小撞擊感度越大［25］。這一規(guī)律適用于離子型或分子型晶體同質異晶體及不同壓力和摻雜條件［26］。新近對硝基苯類化合物的DFT 計算表明，感度主要依賴于LUMO 能級，其次與HOMO 能級有關［27］。

分子中最弱鍵的鍵級或雙原子作用能與炸藥的撞擊感度或熱安定性存在平行的遞變關系，肖鶴鳴等用基態(tài)分子中最弱鍵的鍵級或該鍵連接的雙原子相互作用能來判別同系物的撞擊感度［28］。隨后發(fā)現(xiàn)對引發(fā)鍵為C-NO2的多硝基芳香類炸藥，其Mulliken鍵級與實測感度線性相關。但鍵級判定感度適用于共價型化合物，對于離子型化合物不適用，且只限于單鍵引發(fā)斷裂過程。

DFT 計算表明，化學鍵離解能與相應靜電勢之間存在著關聯(lián)［29］。Politzer等進一步證明，C-NO2鍵中點的靜電勢極值與C-NO2鍵的離解能近似成反比關系，因此用C-NO2鍵中點的靜電勢可以表征C-NO2鍵的強度，進而確定相應硝基化合物的感度［30］。對硝基化合物計算后發(fā)現(xiàn)，熱解引發(fā)的反應活化能Ea與實驗撞擊感度存在線性關系：Ea越小，感度越大［31］。即用熱解反應引發(fā)的活化能判斷結構相似炸藥的撞擊感度［32］。熱解反應引發(fā)的活化能動力學判據(jù)更有普遍性，但局限于分子結構和分解機理類似的化合物。

用DFT 方法對硝基芳香族炸藥計算發(fā)現(xiàn)，感度的對數(shù)與炸藥分子中引發(fā)鍵C-NO2的鍵離解能線性相關［33］。硝基化合物最弱鍵的鍵離解能與整個分子鍵離解能平均值之比（D／Ed，D為最弱鍵的鍵離解能，Ed為分子總的鍵離解能／鍵數(shù)目）可作為沖擊感度的指標，該參數(shù)與16種重要爆炸物的沖擊感度有顯著的線性關系（R≥0.95）；D／Ed值包含整個分子能量因素，在預測沖擊感度方面有更好的表現(xiàn)［34］。

Zeman將多硝基化合物的撞擊與電火花感度，爆炸與熱分解的特征及其13C 及15N 的NMR 化學位移相關聯(lián)［35］。此外，將核獨立化學位移（NICS）與總能量的乘積作為描述符，得到一種新的芳香類炸藥感度評估方法［36］。張朝陽等用硝基靜電荷表示基團的吸電子能力，并將硝基電荷作為判斷硝基芳香族化合物撞擊感度的依據(jù)（即硝基電荷法）［37］。之后又發(fā)現(xiàn)分子中取代基的相互作用能（IE）可判斷硝基化合物衍生物的穩(wěn)定性，并成功地預測硝基苯和硝基苯胺類炸藥的撞擊感度［38］。除了用量子化學參數(shù)描述感度外，近年來分子動力學模擬也用于闡明感度與結構的關系。對TATB 的分子動力學模擬表明，當TATB 固體受到機械撞擊時，其層狀π結構不僅通過分子平面的滑移緩沖機械沖擊作用，還可將沖擊時的部分機械能轉化為分子間作用能，從而減緩分子振動和熱點形成［39］，解釋了TATB的鈍感特性。

對單質和復合含能材料27種樣品的實驗研究表明，特性落高能（H50與落錘質量的乘積）的對數(shù)與熱爆發(fā)活化能之間非常顯著地線性相關［40］。這種適合于不同類型炸藥的關系式，充分表明熱爆發(fā)活化能與感度之間的本質關聯(lián)。如何從理論上求得單質和復合炸藥的熱爆發(fā)活化能，建立既有明確物理意義、又同時適用于單質和復合含能材料的感度預測模型，對理論工作者提出了新的挑戰(zhàn)?？傊?，影響感度的因素很多，僅當引爆過程主要受控于化學作用時，與化學活性相關的量子化學參數(shù)（最小鍵級、分解活化能、前沿分子軌道能隙、離去基團電荷等）才對感度起決定作用。這意味著感度預測模型的多樣性及各自的局限性。

4 分子間相互作用及相容性

無論是由單組分還是多組分構成的含能材料，雖然分子組成和結構對其性能起決定作用，但各組分多聚體或混合組分之間的分子間相互作用，對它們的聚集狀態(tài)、（液體）黏度、（固體）堆積方式和密度以及材料的多種性能（相容性、遷移性能等）也產(chǎn)生重要影響［41］。近10 余年來，量子化學方法已用于研究含能材料分子間的相互作用及其對含能材料性能的影響。對FOX-7二聚體和晶體的DFT 研究發(fā)現(xiàn)，最穩(wěn)定二聚體的構型與晶體中分子堆積方式相似。說明結合能對該晶體中分子的排列方式起決定作用，還預示可按晶體結構較方便地找到某些分子的二聚或多聚體穩(wěn)定構型［42］。在室溫下由FOX-7單體生成最穩(wěn)定二聚體的ΔG＜0，即該二聚體可自發(fā)生成。在常溫下結合能較弱的含能材料通常自發(fā)形成晶體，在很大程度上歸因于分子間相互作用的協(xié)同效應［41］。

分子間相互作用的強弱從本質上決定了多組分體系相容性的大小。對于由分子間作用力結合而成的兩組分體系，當A…A + B…B=2A…B 的ΔG＜0，則A 與B 自發(fā)混合，即完全相容。通常A…B結合能大有利于兩者之間的相容。在實際體系中，A 或B在混合前后均盡可能多地與鄰近分子產(chǎn)生相互作用，并且A 與B分子的大小和形狀各不相同。因此，直接由二聚體結合能值判斷組分間的相容性并不具有普遍性。更為可靠的相容性判斷方法是通過分子動力學模擬，得到溶解度參數(shù)和相互作用參數(shù)等。

從熱力學角度來看，表征體系的相容性可用混合熱ΔHm、混合熵ΔSm、溶解度參數(shù)δ及相互作用參數(shù)。其中溶解度參數(shù)較為簡便，可根據(jù)溶解度參數(shù)差值（Δδ）預測高分子混合物之間的相容性。對于高分子體系，若分子間沒有強極性基團或氫鍵作用，兩種材料的Δδ只要滿足cm-3）1／2兩者就相容［43］。應用分子動力學（MD）和介觀動力學（Meso-Dynamics）對固體推進劑中端羥基聚丁二烯（HTPB）與增塑劑癸二酸二辛酯（DOS）、硝化甘油（NG）的相容性模擬，得到等密度圖、自由能密度和有序度參數(shù)等即可判斷共混體系的相容性；MD 和介觀模擬結果均表明，HTPB／DOS屬于相容體系，DOS 在整個區(qū)域近似均勻分布，與HTPB 相容；而HTPB／NG 屬于不相容體系，與實驗結果一致［43］。

5 固體與界面性質

大多數(shù)含能材料在常溫下是固相。與氣相計算模型相比，固體模型包含鄰近分子的相互影響。用固體模型可研究爆燃傳播、晶格缺陷、機械能與熱能轉化、熱點形成、振動態(tài)激發(fā)與鍵斷裂、晶體能、彈性系數(shù)和膨脹系數(shù)等固體性質。堿金屬疊氮鹽晶體的計算結果表明，其前沿軌道均由疊氮根端位氮的原子軌道組成。由于堿金屬離子對前沿空軌道的貢獻極小，不利于電子從疊氮根向金屬離子躍遷［44］。與此相反，疊氮化銀和疊氮化鉛晶體前沿占有的軌道主要由疊氮根端位氮的原子軌道組成，而前沿空軌道則主要由金屬離子的原子軌道組成，并且?guī)遁^小，這均有助于氮原子上電子直接向金屬離子躍遷［45］。這可解釋堿金屬與過渡金屬疊氮鹽感度的差別。

對于固體計算，人們通常選用計算效率較高的DFT 方法。但DFT 法預測的色散力小于實測值，對分子間作用力較弱的分子晶體，計算結果存在較大誤差。因此，理論化學家致力于對傳統(tǒng)DFT 法進行改進。新近發(fā)展的DFT-D 方法能修正遠程散色作用。用DFT-D 計算10 種典型含能材料分子晶體，其晶胞參數(shù)與實驗結果的相對誤差約2%或更少［46］。在常壓下以及高壓下的晶胞參數(shù)和體積隨壓力變化值與實驗一致，其晶胞參數(shù)的最大誤差為3.67%。不僅能較好地再現(xiàn)分子晶胞結構，還能準確預測晶胞中分子的相對位置和取向，為含能材料分子晶體計算提供更好的結果。

含能材料大多是混合體系，界面現(xiàn)象普遍存在。用DFT 共軛梯度近似（DFT-GGA）研究硝基甲烷和1，1-二胺基-2，2-二硝基乙烯在α-Al2O3（0001）表面上的行為［47］，發(fā)現(xiàn)物理吸附的最穩(wěn)定構型為分子面平行于Al2O3（0001）面。硝基甲烷分解的最低能壘（58.5kJ／mol）為H 原子消去；其C-N斷裂能壘為158.8kJ／mol，均小于氣相分解時的能壘。FOX-7在α-Al2O3（0001）表面上不發(fā)生分解。表明Al2O3表面氧化層雖有純化作用，但對炸藥分解有一定的催化能力［47］。用DFT-GGA 法研究硝胺分子在Al（111）表面的吸附表明［48］：既有物理吸附又有化學吸附。物理吸附時，硝基氧原子指向鋁表面。化學吸附時，氧和氮原子與表面鋁原子形成化學鍵，N-O 和N-N 鍵斷裂，形成強的Al-O 和Al-N 鍵，鋁表面可被硝基和氨基分解出的N 和O原子迅速氧化。DFT-GGA 法研究TNT 在鋁（111）表面吸附結果表明，氧和鋁原子間強烈的吸引作用使TNT 中N-O 鍵分解。TNT 鄰位硝基的N-O 鍵比對位硝基的N-O 鍵易斷裂。除了硝基的N-O 鍵斷裂外，其他鍵未斷裂，這有別于TNT 直接分解時C-NO2鍵優(yōu)先斷裂?？梢娕cAl表面的作用改變了反應歷程。含能材料在Al表面的分解過程，實際上是含鋁炸藥爆燃后期，Al表面氧化層破裂后發(fā)生的反應。因此，綜合考慮金屬氧化層及其厚度、氧化層高溫破裂過程對含能材料吸附和分解的影響尤為重要。

固體和界面的一些性質，特別是涉及動態(tài)過程時，需要借助分子動力學模擬。用第一性原理分子動力學模擬RDX 在Al（111）表面的吸附和分解的結果表明，氧與鋁的強吸引作用使RDX 的N-O和N-N 鍵斷裂并產(chǎn)生NO；此外，RDX 環(huán)也發(fā)生斷裂［49］。分子動力學模擬HMX 與石墨界面的彎折和滑移行為，得到表面結構和能量變化過程，表面弛豫和摩擦均對熱點形成有貢獻［50］。采用MD模擬鍵合劑對丁羥推進劑中端羥基聚丁二烯（HTPB）與Al／Al2O3之間界面的吸附能與力學性能的影響時發(fā)現(xiàn)［51］：鍵合劑在Al2O3晶面的吸附能高于HTPB在Al2O3晶面的吸附能，而在Al晶面的規(guī)律并不明顯。Al／HTPB吸附界面加入鍵合劑TEA后，隨其濃度的增加，體系彈性模量增大。對DOS、IPDI和HTPB組成的丁羥推進劑黏接體系的運動過程進行模擬［52］，求得增塑劑DOS 在黏合劑體系（HTPB+IPDI）中的擴散系數(shù)，擴散系數(shù)隨環(huán)境溫度的升高而增大，隨增塑劑含量的增大而減小，模擬結果與實驗結果吻合，揭示了微觀遷移現(xiàn)象。由于微觀或介觀模型的建立及模擬條件的選取常存在一定隨意性，影響模擬結果的再現(xiàn)性，故特別強調模型應包含真實體系主要特征并且有統(tǒng)計意義。

6 極端條件下的性質

在實驗技術難以實現(xiàn)的極端高壓和高溫下，用量子化學計算含能材料高壓性質和MD 模擬其化學演變過程顯得尤為重要。但經(jīng)典MD不能對化學反應體系進行模擬，近10年發(fā)展的反應力場（ReaxFF）方法［53］和從頭算CPMD方法為模擬含能材料高溫高壓反應過程提供了新的技術途徑。ReaxFF反應力場模型依據(jù)鍵距與鍵級、鍵級與鍵能的關系，描述化學鍵斷裂時體系能量變化，在原來的經(jīng)典力場中加入鍵級參數(shù)及其修正項，比量子化學計算速度快萬倍［53］，可以模擬含數(shù)百萬原子體系的反應過程［54］。

硝基甲烷在均勻受壓和3個軸向受壓時的量子化學計算結果表明［55］，在靜水壓下，當體積壓縮至原來的50%時，壓力上升至50GPa，ΔEHOMO-LUMO減小了約0.6eV。最高占有軌道（HOMO）和最低空軌道（LUMO）同時增加，其差值卻幾乎隨體積的變化單調下降，帶隙減小從理論上可以認為感度在增加。在y方向上，帶隙最多可減小1.4eV，比靜水壓時帶隙下降值大，并且?guī)冻史菃握{變化。x方向的應變引起的能隙變化與y方向的類似，只不過當V／V0接近50%時y方向帶隙會引起陡降。當V／V0為65%～80%而z方向的能隙急劇下降。z方向壓縮時，分子相互靠近時會導致一個分子的甲基與其鄰近分子的硝基靠近，定域在硝基上HOMO和LUMO 的電子密度比其他類型的壓縮形變更大，能隙下降得更快。高壓所產(chǎn)生的y軸向張力，導致單胞內每個硝基甲烷分子均有1個C-H 鍵產(chǎn)生急劇伸展。隨著y軸應變增加至50%，會導致這些鍵被進一步拉伸并引起質子消去。

HMX 和TATB 高溫分解的分子動力學模擬表明，TATB經(jīng)快速分解（30ps）產(chǎn)生大量多芳環(huán)為主的碳簇合物（總量的15%～30%），HMX 則產(chǎn)生小分子產(chǎn)物。HMX 在較低溫度下易于分解，TATB分解速率比HMX 小一個數(shù)量級［54］。這可解釋兩種炸藥分解過程實驗現(xiàn)象的差別。分子模擬還揭示高溫分解時各物種濃度變化。

對液體和固體硝基甲烷在2 000～3 000K、密度1.97g／cm3、時間200ps內的高溫過程進行模擬［56］，在3 000K 時分解起始步驟為分子間質子遷移，生成CH3NOOH 和CH2NO2。在較低的溫度（2 500～2 000K）第一步反應則為異構化生成CH3ONO 的過程。同時在反應剛開始極短時間內可伴隨分子內質子遷移生成CH2NOOH。作為放熱過程的標志，H2O 的形成則是在這些過程之后。2000K 時液態(tài)硝基甲烷異構化生成亞硝酸酯，與從頭算結果相一致。

采用MD 研究含能材料的分解過程，一般用ReaxFF力場。目前的ReaxFF參數(shù)通過擬合DFT計算結果得到。而DFT 方法低估了分子間的London色散力，導致晶體平衡體積偏高約10%～15%。最近通過對長程色散作用進行低梯度校正，得到ReaxFF-lg力場［57］。由室溫下RDX 分子質心距離徑向分布函數(shù)可知，ReaxFF-lg使分子模擬結果大為改進。預期這類修正的反應力場將對高溫高壓條件的含能材料研究發(fā)揮更有效的作用。

7 結束語

（1）隨著量子化學和分子動力學模擬方法在含能材料領域的廣泛應用，研究對象逐漸從分子結構簡單體系擴展到復雜體系，并從單相擴展到多相、從體相擴展到界面。

（2）分子動力學模擬從平衡態(tài)模擬發(fā)展到施加外場的非平衡態(tài)模擬。尤其是介觀層次的模擬不僅揭示納米材料的特殊性質，還構筑了微觀結構與宏觀性質之間的橋梁。但介觀體系涉及量子效應，如何改進量子化學方法，以實現(xiàn)介觀層次的量子化學研究，對理論工作者提出了挑戰(zhàn)。

（3）可以預測，量子化學和分子模擬對開發(fā)新型含能材料將發(fā)揮明顯的、不可替代的作用。但是理論模擬不能代替實驗，可靠的實驗數(shù)據(jù)是檢驗模擬模型正確性的唯一依據(jù)。模擬和實驗是相輔相成、缺一不可的研究手段，兩者的結合將極大地促進含能材料的發(fā)展。

致謝：本綜述的撰寫得益于燃燒與爆炸技術重點實驗室趙鳳起研究員的鼓勵、支持、并提出寶貴的修改意見，在此深表感謝！

［1］Repasky M P，Chandrasekhar J，Jorgensen W L.PDDG／PM3and PDDG／MNDO：Improved semiempirical methods［J］.J Comput Chem，2002，23：1601-1622.

［2］Dorsett H，White A.Overview of molecular modelling and ab initio molecular orbital methods suitable for use with energetic materials［R］.Salisbury：DSTO Aeronautical and Maritime Research Laboratory，2000.

［3］Keshavarz M H.Improved prediction of heats of sublimation of energetic compounds using their molecular structure［J］.J Hazard Mater，2010，177（1-3）：648-659.

［4］Gharagheizi F，Sattari M，Tirandazi B.Prediction of crystal lattice energy using enthalpy of sublimation：A group contribution-based model［J］.Ind Eng Chem Res，2011，50：2482-2486.

［5］Kamlet M J，Jacobs S J.Chemistry of detonations.I.A simple method for calculating detonation properties of C-H-N-O explosives［J］.J Chem Phys，1968，48：23-35.

［6］邱玲.氮雜環(huán)硝胺類高能量密度材料（HEDM）的分子設計［D］.南京：南京理工大學，2007.

QIU Ling.Molecular design of high energy density materials（HEDM）-cyclic nitramines［D］.Nanjing：Nanjing University of Science and Technology，2007.

［7］Ju X H，Li Y M，Xiao H M.Theoretical studies on the heats of formation and the interactions among the difluoroamino groups in polydifluoroaminocubanes［J］.J Phys Chem A，2005，109：934-938.

［8］Fan X W，Ju X H，Xia Q Y，et al.Strain energies of cubane derivatives with different substituent groups［J］.J Hazard Mater，2008，151（1）：255-260.

［9］Fan X W，Ju X H.Theoretical studies on four-membered ring compounds with NF2，ONO2，N3，and NO2groups［J］.J Comput Chem，2008，29：505-513.

［10］馮麗娜，張建國，張同來，等.氨基四唑化合物異構和分解反應的研究進展［J］.含能材料，2009，17：113-118.

FENG Li-na，ZHANG Jian-guo，ZHANG Tong-lai，et al.Progress in the tautomerism and decomposition of amino-tetrazoles［J］.Chin J Energetic Mater，2009，17：113-118.

［11］霍冀川，吳瑞榮，舒遠杰，等.氮雜環(huán)類含能材料熱分解研究進展［J］.爆破，2007，24：21-25.

HUO Ji-chuan，WU Rui-rong，SHU Yuan-jie，et al.Investigative development of thermal decomposition of zaohetercycline energetic materials［J］.Blasting，2007，24：21-25.

［12］王新鋒，舒遠杰.新型高能炸藥熱分解研究進展［J］.化學研究與應用，2004，16（3）：305-308.

WANG Xin-feng，SHU Yuan-jie.Studies on thermal decomposition of new high explosives［J］.Chem Res Appl，2004，16（3）：305-308.

［13］Alavi S，Reilly L M，Thompson D L.Theoretical predictions of the decomposition mechanism of 1，3，3-trinitroazetidine（TNAZ）［J］.J Chem Phys，2003，119（16）：8297-8304.

［14］Zhao Q，Zhang S，Li Q S.A direct ab initio dynamics study of the initial decomposition steps of gas phase 1，3，3-trinitroazetidine［J］.Chem Phys Lett，2005，412：317-321.

［15］Kohno Y，Takahashi O，Saito K.Theoretical study of initial decomposition process of NTO dimer［J］.Phys Chem Chem Phys，2001，3：2742-2746.

［16］Wang Y M，Chen C，Lin S T.Theoretical studies of the NTO unimolecular decomposition［J］.J Mol Struct（THEOCHEM），1999，460（1）：79-102.

［17］Kiselev V G，Gritsan N P.Theoretical study of the 5-aminotetrazole thermal decomposition［J］.J Phys Chem A，2009，113（15）：3677-3684.

［18］Guo Y Q，Greenfield M，Bernstein E R.On the excited electronic state dissociation of nitramine energetic materials and model systems［J］.J Chem Phys，2007，127：154301.

［19］Bhattacharya A，Bernstein E R.Nonadiabatic decomposition of gas-phase RDX through conical intersections：An ONIOM-CASSCF study［J］.J Phys Chem A，2011，115：4135-4147.

［20］Bhattacharya A，Guo Y Q，Bernstein E R.Nonadiabatic reaction of energetic molecules［J］.Acc Chem Res，2010，43（12）：1476-1485.

［21］Zhang L，Zybin S V，van Duin A C T，et al.Shock induced decomposition and sensitivity of energetic materials by ReaxFF molecular dynamics.Shock Compression of Condensed Matter［M］.New York：American Institute of Physics，2006：585.

［22］Manaa M.R，Reed E J，F(xiàn)ried L E，Goldman N.Nitrogen-rich heterocycles as reactivity retardants in shocked insensitive explosives［J］.J Am Chem Soc，2009，131：5483-5487.

［23］Kamlet M J，Adolph H G.The relationship of impact sensitivity with structure of organic high？explosives II.Polynitroaromatic explosives［J］.Prop，Expl，Pyro，1979，4（2）：30-34.

［24］Delpuech A，Cherville J.Relation entre la structure electronique et la sensibilitéau choc des explosifs secondaires nitrés.Critère moléculaire de sensibilitéII.Cas des esters nitriques［J］.Prop，Expl，Pyro，1979，4（6）：121-128.

［25］肖鶴鳴.金屬疊氮化物的能帶和電子結構［M］.北京：科學出版社，1996.

［26］朱衛(wèi)華，張效文，肖鶴鳴.高能晶體撞擊感度理論研究：第一性原理帶隙（ΔEg）判據(jù)［J］，含能材料，2010，18（4）：431-434.

ZHU Wei-hua，ZHANG Xiao-wen，XIAO He-ming.Theoretical studies of impact sensitivity of energetic crystals：First-principles band gap（ΔEg）criterion［J］.Chin J Energetic Mater，2010，18（4）：431-434.

［27］Türker L.Tunneling effect and impact sensitivity of certain explosives［J］.J Hazard Mater，2009，169：819-823.

［28］肖鶴鳴.硝基化合物的分子軌道理論［M］.北京：國防工業(yè)出版社，1993.

［29］March N H.The role of the bond midpoint electron density in homonuclear molecular binding［J］.Int J Quantum Chem，1994，52（1）：247-265.

［30］Politzer P，Murray J S.Relationships between dissociation energies and electrostatic potentials of C-NO2bonds：Applications to impact sensitivities［J］.J Mol Struct（THEOCHEM），1996，376（1-3）：419-424.

［31］Fan J F，Gu Z M，Xiao H M，et al.Theoretical study on pyrolysis and sensitivity of energetic compounds：（4）Nitro derivatives of phenols［J］.J Phys Org Chem，1998，11（3）：177-184.

［32］Chen Z X，Xiao H M.Impact sensitivity and activation energy of pyrolysis for tetrazole compounds［J］.Int J Quantum Chem，2000，79：350-357.

［33］Rice B M，Sahu S，Owens F J.Density functional calculations of bond dissociation energies for NO2scission in some nitroaromatic molecules［J］.J Mol Struct（THEOCHEM），2002，583：69-72.

［34］Mathieu D.Theoretical shock sensitivity index for explosives［J］.J Phys Chem A，2012，116（7）：1794-1800.［35］Zeman S.New Aspects of impact reactivity of polynitro compounds，Part III.Impact sensitivity as a function of the intermolecular interactions［J］.Prop，Expl，Pyro，2003，28（6）：301-307.

［36］杜軍良，舒遠杰，周陽.一種表征芳香族炸藥撞擊感度的簡單方法［J］.含能材料，2008，16（6）：766-766.

DU Jun-liang，SHU Yuan-jie，ZHOU Yang.A simple method for predicting impact sensitivity of aromatic explosives［J］.Chin J Energetic Mater，2008，16（6）：766-766.

［37］Zhang C Y，Shu Y J，Huang Y G，et al.Investigation of correlation between impact sensitivities and nitro group charges in nitro compounds［J］.J Phys Chem B，2005，109：8978-8982.

［38］Zhang C Y.Investigation on the correlation between the interaction energies of all substituted groups and the molecular stabilities of nitro compounds［J］.J Phys Chem A，2006，110：14029-14035.

［39］Zhang C Y，Wang X C，Huang H.π-Stacked interac-tions in explosive crystals：Buffers against external mechanical stimuli［J］.J Am Chem Soc，2008，130：8359-8365.

［40］劉子如，岳璞，任曉寧，等.熱爆發(fā)活化能研究［J］.火炸藥學報，2011，34（6）：58-63.

LIU Zi-ru，YUE Pu，REN Xiao-ning，et al.Investigation on activation energy of heat explosion［J］.Chinese Journal of Explosives and Propellants，2011，34（6）：58-63.

［41］肖鶴鳴，居學海.高能體系中的分子間相互作用［M］.北京：科學出版社，2004.

［42］Ju X H，Xiao H M，Xia Q Y.A density functional theory investigation of 1，1-diamino-2，2-dinitroethylene dimers and crystal［J］.J Chem Phys，2003，119（19）：10247-10255.

［43］趙貴哲，馮益柏，付一政，等.端羥基聚丁二烯／增塑劑共混物相容性的分子動力學模擬和介觀模擬［J］.化學學報，2009，67（19）：2233-2238.

ZHAO Gui-zhe，F(xiàn)ENG Yi-bai，F(xiàn)U Yi-zheng，et al.Molecular dynamic simulations and mesoscopic dynamic simulations on the compatibility of HTPB／plasticizer blends［J］.Acta Chim Sin，2009，67（19）：2233-2238.

［44］Ju X H，Xiao H M，Guang F J.Periodic DFT approaches to crystalline alkali metal azides［J］.Chin Sci Bull，2002，47：1180-1183.

［45］居學海，姬廣富，邱玲.Cu+、Ag+疊氮鹽晶體的周期性ab initio 計 算［J］.高 等 學 校 化 學 學 報，2005，26

（11）：2125-2127.JU Xue-hai，JI Guang-fu，QIU Ling，et al.Periodic ab initio calculations on Ag+and Cu+azides［J］.Chem J Chin Univ，2005，26（11）：2125-2127.

［46］Sorescu D C，Rice B M.Theoretical predictions of energetic molecular crystals at ambient and hydrostatic compression conditions using dispersion corrections to conventional density functionals（DFT-D）［J］.J Phys Chem C，2010，114（14）：6734-6748.

［47］Sorescu D C，Boatz J A，Thompson D L.First-principles calculations of the adsorption of nitromethane and 1，1-diamino-2，2-dinitroethylene （FOX-7）molecules on theα-Al2O3（0001）surface［J］.J Phys Chem B，2005，109（4）：1451-1463.

［48］Zhou S Q，Zhao F Q，Ju X H，et al.A density functional theory study of adsorption and decomposition of nitroamine molecules on the Al（111）surface［J］.J Phys Chem C，2010，114（20）：9390-9397.

［49］Umezawa N，Kalia R K，Nakano A，et al.1，3，5-Trinitro-1，3，5-triazine decomposition and chemisorption on Al（111）surface：First-principles molecular dynamics study［J］.J Chem Phys，2007，126：234702.

［50］Long Y，Liu Y G，Nie F D，et al.Theoretical study of breaking and slipping processes for HMX／graphite interface［J］.Appl Surf Sci，2012，258（7）：2384-2392.

［51］焦東明，楊月誠，強洪夫，等.鍵合劑對HTPB 與Al／Al2O3之間界面作用的分子模擬［J］.火炸藥學報，2009，32（4）：60-763

JIAO Dong-ming，YANG Yue-cheng，QIANG Hongfu，et al.Molecular simulation of effect of bonding agents on interface interaction for HTPB and Al／Al2O3［J］.Chinese Journal of Explosives and Propellants，2009，32（4）：60-763.

［52］李紅霞，強洪夫，武文明.丁羥推進劑黏結體系中增塑劑遷移的分子模擬［J］.火炸藥學報，2008，31（5）：74-78.

LI Hong-xia，QIANG Hong-fu，WU Wen-ming.Molecular simulation on plasticizer migration in the binding system of HTPB propellant［J］.Chinese Journal of Explosives and Propellants，2008，31（5）：74-78.

［53］van Duin A C T，Dasgupta S，Lorant F，et al.ReaxFF：A reactive force field for hydrocarbons［J］.J Phys Chem A，2001，105：9396-9409.

［54］Zhang L，Zybin S V，van Duin A C T，et al.Carbon cluster formation during thermal decomposition of octahydro-1，3，5，7-tetranitro-1，3，5，7-tetrazocine and 1，3，5-Triamino-2，4，6-trinitrobenzene high explosives from ReaxFF reactive molecular dynamics simulations［J］.J Phys Chem A，2009，113：10619-10640.

［55］Manaa M R，F(xiàn)ried L E，Reed E J.Explosive chemistry：Simulating the chemistry of energetic materials at extreme conditions［J］.J Comput Aid Mol Des，2003，10（2）：75-97.

［56］Han S P，van Duin A C T，Goddard III W A，et al.Thermal decomposition of condensed-Phase nitromethane from molecular dynamics from ReaxFF reactive dynamics［J］.J Phys Chem B，2011，115：6534-6540.

［57］Liu L C，Liu Y，Zybin S V，et al.Correction of the ReaxFF reactive force field for London dispersion，with applications to the equations of state for energetic materials［J］.J Phys Chem A，2011，115：11016-11022.