王麗莉，熊 鷹，謝煒宇，牛亮亮，張朝陽

（1.中國工程物理研究院計(jì)算機(jī)應(yīng)用研究所，四川 綿陽 621999；2.中國工程物理研究院化工材料研究所，四川 綿陽 621999）

1 引言

含能材料密度的理論預(yù)測方法可分成以下三大類：一是基于分子體積計(jì)算密度，由于分子體積沒有明確的物理定義，因此存在多種計(jì)算方法，如等電子密度面（Isosurface of Electron Density，IED）方法通常將分子摩爾體積定義為0.00l a.u.電子密度等值面輪廓范圍內(nèi)的體積［1］，分子表面靜電勢（Molecular Surface Electrostatic Potentials，MESP）校正法是在密度求解公式中加入通過實(shí)驗(yàn)數(shù)據(jù)擬合的系數(shù)表征分子間相互作用，基團(tuán)加和（Group Additivity Methods，GAM）法是由原子和官能團(tuán)體積的貢獻(xiàn)估算分子體積等，其中計(jì)算電子密度的方法多數(shù)采用密度泛函理論（Density Functional Theory，DFT）、二階微擾（Second-order Moller-Plesset Perturbation，MP2）理論、耦合簇（Coupled Cluster with Singles and Doubles，CCSD）理論等。二是基于晶體體積計(jì)算密度，晶體結(jié)構(gòu)預(yù)測（Crystal Structure Prediction，CSP）結(jié)合了能量和結(jié)構(gòu)的精確計(jì)算方法與全局結(jié)構(gòu)搜索方法，其中結(jié)構(gòu)與能量的計(jì)算方法包括密度泛函理論的色散校正（Dispersion Correction Density Functional Theory，DFT-D），分子力場與分子動力學(xué)方法（Molecular Dynamics，MD），半經(jīng)驗(yàn)分子軌道方法（Semi-empirical Molecular Orbital，MO）等；結(jié)構(gòu)搜索的代表性算法有模擬退火法（Simulated Annealing，SA）、隨機(jī)結(jié)構(gòu)搜索法（Ab Initio Random Structure Searching，AIRSS）、進(jìn)化算法（Evolutionary Approach，EA）和粒子群優(yōu)化算法（Particle Swarm Optimization，PSO）等。三是基于數(shù)學(xué)和統(tǒng)計(jì)方法估算含能材料密度，如通過多元線性回歸（Multiple Linear Regressin，MLR）、人工神經(jīng)網(wǎng)絡(luò)（Artificial Neural Networks，ANN）和改進(jìn)的遺傳算法（Genetic Algorithm，GA）等建立快速穩(wěn)定的定量構(gòu)效關(guān)系（Quantitative Structure Property Relationship，QSPR）模型，或根據(jù)經(jīng)驗(yàn)公式由元素組成和官能團(tuán)數(shù)目估算晶體密度等。隨著云計(jì)算的普及和大數(shù)據(jù)時(shí)代的到來，機(jī)器學(xué)習(xí)（Machine Learning，ML）及人工智能（Artificial Intelligence，AI）技術(shù)將推動含能晶體結(jié)構(gòu)、性質(zhì)以及感度之間關(guān)系的研究，并為開發(fā)新型鈍感高能炸藥提供依據(jù)。

2 基于分子體積計(jì)算密度

2.1 等電子密度面法

等電子密度面法計(jì)算密度的主要步驟如下：優(yōu)化含能分子確保其處于能量最低點(diǎn)，然后計(jì)算振動光譜和相應(yīng)熱力學(xué)信息評估結(jié)構(gòu)穩(wěn)定性［2］，用Monte-Carlo方法計(jì)算每個(gè)分子的平均摩爾體積V（0.001），理論密度用分子摩爾質(zhì)量M與V（0.001）之比求得［3］。該方法對分子形狀和體積的描述依賴電子密度值的選擇，對于處于凝聚態(tài)的分子，考慮到分子間壓積會使范德華表面有一定穿透，也有用0.002 a.u.或0.0015 a.u.的電子密度等值面作為范德華表面［4-6］。

密度泛函理論的B3LYP 方法和6-31G**基組是等電子密度面法計(jì)算分子晶體密度的主要方法［7］，準(zhǔn)確性優(yōu)于半經(jīng)驗(yàn)MO 方法［8］，經(jīng)濟(jì)性優(yōu)于MP2 和CCSD（T）等高水平電子密度方法［9］。對CHNO 型分子晶體，該方法計(jì)算中性分子的理論密度與實(shí)驗(yàn)值的平均誤差和均方根誤差分別為1%和3.7%，分子中含氮量超過50%或含氮環(huán)的誤差為-1.8%和3.4%，離子晶體的誤差為-5.1%和6.6%，通過對體積公式進(jìn)行修正，離子晶體平均誤差和均方根誤差減小至-1.1%和4.8%［10］。

等電子密度面法沒有考慮分子間空隙及各種復(fù)雜的相互作用，用于密度這種宏觀性質(zhì)計(jì)算，原理受質(zhì)疑。由于該方法對極性分子的計(jì)算誤差偏大，后期根據(jù)經(jīng)驗(yàn)進(jìn)行參數(shù)的計(jì)算擬合，對不同的分子結(jié)構(gòu)進(jìn)行體積修正，精確度有所提高。

2.2 分子表面靜電勢校正法

P.Politzer 等［11-12］提出修正分子表面靜電勢參數(shù)法，在密度求解公式中分別考慮了與分子表面正、負(fù)電荷分布情況相關(guān)的分子表面靜電勢平衡系數(shù)v和總方差計(jì)算公式見式（1）：

式中，M為分子摩爾質(zhì)量，g·mol-1；α、β、γ是衰減系數(shù)，由計(jì)算擬合得到。該方法預(yù)測密度與實(shí)驗(yàn)值誤差一般小于0.05 g·cm-3，被廣泛應(yīng)用于CHNO 型化合物密度的預(yù)測［13］。如H.Singh［14］、P.Ma［15］和Y.J.Shu等［16］用Politzer 方法計(jì)算了2，4，6-三硝基甲苯（TNT）、環(huán)三亞甲基三硝胺（RDX）、環(huán)四亞甲基四硝胺（HMX）等多種含能分子的結(jié)構(gòu)和密度，都與實(shí)驗(yàn)數(shù)據(jù)吻合較好。2013 年，B.M.Rice［17］用Politzer 方法對前期工作進(jìn)行改進(jìn)［10］，并將計(jì)算結(jié)果分組與實(shí)驗(yàn)數(shù)據(jù)比較，其中38 個(gè)中性分子晶體理論密度的平均絕對誤差從0.05 g·cm-3減小至0.035 g·cm-3，48 個(gè)離子化合物理論密度的平均絕對誤差從0.088 g·cm-3減小至0.045 g·cm-3。張君君等［18］用Politzer 方法計(jì)算出N-1，4，6-三硝基六氫咪唑（TNINA）和［4，5-d］咪唑-2（1H）-亞硝胺（DNINA）晶體密度分別為1.91 g·cm-3和1.83 g·cm-3，與實(shí)驗(yàn)密度1.89 g·cm-3和1.82 g·cm-3吻合較好。C.K.Kim［19］將密度計(jì)算公式修正為如下形式：

式中，AREA 是范德華分子表面積，nm2；Π 是平均偏差，是靜電勢為分子表面上平均電勢，kJ·mol-1；該方法計(jì)算密度的平均絕對誤差是0.039 g·cm-3。A.Nirwan 等［20］將221 種炸藥分成7 類，對比了Lee、Rice、Kim、Politzer 等MESP 方法計(jì)算密度的有效性，其中Rice 和Politzer 修正更接近實(shí)驗(yàn)值，多數(shù)誤差小于3%，硝酸脂和包含籠狀或剛性環(huán)結(jié)構(gòu)的誤差較大。

分子表面靜電勢校正方法是目前計(jì)算含能材料密度的主要方法，計(jì)算精度優(yōu)于等電子密度面法和基團(tuán)加和法，其可靠性己經(jīng)被證明。

2.3 基團(tuán)加和法

基團(tuán)加和法用組成分子的原子和官能團(tuán)的標(biāo)準(zhǔn)體積之和估算分子體積，再用分子摩爾質(zhì)量求密度。H.L.Ammon 等［21］從11555 個(gè)晶體結(jié)構(gòu)數(shù)據(jù)中提取了78 種原子和官能團(tuán)的體積，建立標(biāo)準(zhǔn)體積數(shù)據(jù)庫，用于計(jì)算CHNOF 類化合物的晶體密度，后來又?jǐn)U充和升級該數(shù)據(jù)庫［22-23］，用于更多有機(jī)物的密度預(yù)測。D.Mathieu［24］以Ammon 數(shù)據(jù)庫為基礎(chǔ)，用基團(tuán)加和法預(yù)測密度與實(shí)驗(yàn)數(shù)據(jù)的平均誤差接近2%，只有部分密堆積結(jié)構(gòu)的估算密度誤差達(dá)到6%～8%。S.Beaucamp［25-26］用帶電基團(tuán)體積加和法估算離子鹽和水合物的晶體密度，平均誤差小于2.5%，通過修正氫鍵和環(huán)狀結(jié)構(gòu)等對分子體積的貢獻(xiàn)［27］，42880 個(gè)晶體估算密度的平均誤差為2%。

基團(tuán)加和法中體積的貢獻(xiàn)值是從大量實(shí)驗(yàn)數(shù)據(jù)中總結(jié)得出，缺乏實(shí)驗(yàn)數(shù)據(jù)難以獲得準(zhǔn)確的經(jīng)驗(yàn)參數(shù)。不僅如此，該方法不考慮分子構(gòu)型和分子間相互作用，無法區(qū)分同分異構(gòu)體的密度差異，忽略晶型及溫度對密度的影響，對非電中性含能化合物，含硅、硼等元素的化合物，或部分密堆積結(jié)構(gòu)含能晶體的計(jì)算偏差較大，估算范圍有一定限制。

3 基于晶體結(jié)構(gòu)計(jì)算密度

3.1 量化方法優(yōu)化晶體結(jié)構(gòu)

從理論上預(yù)測含能晶體的結(jié)構(gòu)，可以得到密度。晶體結(jié)構(gòu)的求解通常需要高精度量化方法進(jìn)行結(jié)構(gòu)優(yōu)化來實(shí)現(xiàn)局域最小值精修和能量穩(wěn)定性計(jì)算。

3.1.1 密度泛函理論的色散校正

密度泛函理論對分子間弱相互作用的描述不完全準(zhǔn)確，為了能夠精確的計(jì)算含能分子晶體的物性，需要加入分子間的色散校正。非局域性色散校正、London 色散校正、球原子模型等是對DFT 的有效修正。其中DFT-D 方法在量化程序中得到廣泛支持，通常用零阻尼函數(shù)（式（5））來調(diào)節(jié)色散校正在近程、中程距離時(shí)的行為，避免近距離時(shí)校正能（式（4））較大導(dǎo)致的重復(fù)問題［28］。

其中RAB代表AB 原子間的距離，sn是截?cái)喟霃?，nm；C是原子間色散校正系數(shù)；sn和sr，n是刻度因子；γ是預(yù)設(shè)常數(shù)。隨原子間距離減小，f逐漸為0，故稱零阻尼。

經(jīng)過DFT-D 校正后，靜電主導(dǎo)的弱相互作用，如氫鍵、鹵鍵等問題，計(jì)算精度有較大改進(jìn)，已用于計(jì)算5，5′-聯(lián)四唑-1，1′-二氧二羥銨鹽（TKX-50）、1，3，5-三氨基-2，4，6-三硝基苯（TATB）、1，1-二氨基-2，2-二硝基乙烯（FOX-7）、2，6-二氨基-3，5-二硝基吡嗪-1-氧（LLM-105）、β-HMX 和六硝基六氮雜異戊茲烷（CL-20）等晶體密度［29-30］，晶胞參數(shù)計(jì)算結(jié)果與實(shí)驗(yàn)的誤差通常在2%以內(nèi)［31］。L.Zhang等［32］采用PBE交換關(guān)聯(lián)勢計(jì)算HMX 的α、β、δ和γ相的晶體結(jié)構(gòu)，體積誤差分別為-2.36%，-1.00%，+0.12%和-0.3%，采用vdW-D3方法計(jì)算了53 種炸藥晶體結(jié)構(gòu)，體積誤差＜±4%［33］。H.Lin［34］、H.C.S.Chan［35］、P.Panini［36］也采用DFT-D對含能晶體或共晶進(jìn)行結(jié)構(gòu)和密度預(yù)測。

DFT-D 是計(jì)算含能材料晶體密度最適合的DFT方法，計(jì)算精度由泛函和色散校正模型的優(yōu)良性、兩者對色散作用描述的互補(bǔ)性兩個(gè)因素決定，誤差明顯小于傳統(tǒng)的DFT 方法。DFT-D 的可行性與精度經(jīng)過了系統(tǒng)和全面的驗(yàn)證，但計(jì)算量隨體系的增大呈三次方或二次方增長，比較耗時(shí)。采用對稱適應(yīng)微擾理論（Symmetry Adapted Perturbation Theory，SAPT）把相互作用能直接分解為靜電相互作用、誘導(dǎo)能和色散能、對稱性微擾理論定義的交換能，可以計(jì)算更大規(guī)模的原子模型。R.Podeszwa［37］采用基于DFT 單體描述的SAPT 準(zhǔn)確預(yù)測了RDX 的晶體結(jié)構(gòu)，計(jì)算了FOX-7 的完整勢能面。

3.1.2 分子動力學(xué)結(jié)合第一原理

提高DFT 方法計(jì)算規(guī)模和效率的方法之一是與分子力場和分子動力學(xué)方法結(jié)合。MD 方法可研究單質(zhì)炸藥及其復(fù)合體系的結(jié)構(gòu)、晶體生長機(jī)制［38］、能量、熱穩(wěn)定性、堆積密度、爆轟性能、溶劑對晶體的影響、力學(xué)性能及其溫度效應(yīng)，并與感度關(guān)聯(lián)研究其安全性能，計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)越來越接近［39］。

分子動力學(xué)結(jié)合第一原理計(jì)算在晶體密度求解方面取得了極大的進(jìn)展。P.Srinivasan［40］用B3LYP 方法和6-311G*基組優(yōu)化2，4-二硝基苯甲酸（DNBA）的分子結(jié)構(gòu)，再用MOLPAK（MOLecular PAcKing）建立空間堆積模型，進(jìn)行UMD 力場晶格能最小化［41］，晶體結(jié)構(gòu)預(yù)測結(jié)果與X-ray 衍射實(shí)驗(yàn)結(jié)果的偏差很小，密度的計(jì)算結(jié)果與實(shí)驗(yàn)值誤差為2.98%。D.N.Kaushik［42］用第一原理處理有機(jī)晶體中的單分子和短程二體相互作用，用極化分子力場處理多體和長程相互作用，建立了基于片段的雜化多體相互作用模型，晶胞參數(shù)計(jì)算結(jié)果優(yōu)于色散修正B3LYP-D*和AMOEBA 極化力場方法。J.S.Li［43］在B3LYP 方法和6-31++G**基組計(jì)算分子結(jié)構(gòu)的基礎(chǔ)上，用Dreiding 力場和模擬退火法預(yù)測了2，4，6，8-四硝基-1，3，5，7-四氮雜環(huán)丁烷（TNTAzC）和2，3，5，6，7，8-六硝基-1，4-二氮雜環(huán)丁烷（HNDAzC）的晶體密度分別為2.156 g·cm-3和1.975 g·cm-3。力場選擇對晶體堆積和密度計(jì)算非常關(guān)鍵，Dreiding 力場對84%的不含芳香環(huán)結(jié)構(gòu)的密度預(yù)測值與實(shí)驗(yàn)值的偏差在5%以內(nèi)，PCFF 力場對83%的含芳香環(huán)結(jié)構(gòu)的密度預(yù)測值與實(shí)驗(yàn)值的偏差在4%以內(nèi)［44］，廣義Amber 力場計(jì)算FOX-7、RDX、CL-20、HMX 等的密度比實(shí)驗(yàn)值低5%～8%［45］，CVFF 力場和DMACRYS 力場也用于含芳香環(huán)結(jié)構(gòu)的含能晶體預(yù)測，密度誤差小于7%［46-47］。

分子動力學(xué)的核心在于求解正則運(yùn)動方程，分子內(nèi)和分子間的作用力都包含在分子力場中，因此計(jì)算的準(zhǔn)確性依賴所使用的力場［48-49］。如果缺乏完全適用的力場參數(shù)，無法準(zhǔn)確描述分子間的相互作用，致使晶體結(jié)構(gòu)預(yù)測存在較大偏差，要提高晶體密度計(jì)算的準(zhǔn)確性需自研發(fā)經(jīng)驗(yàn)勢。

3.1.3 半經(jīng)驗(yàn)方法

半經(jīng)驗(yàn)方法通過忽略分子積分計(jì)算或引進(jìn)經(jīng)驗(yàn)參量對Hartree-Fock 方程的近似求解，如電荷自恰密度泛函緊束縛（Self-consistent Charge Density Functional Tight-binding，SCC-DFTB）方法是在DFT 基礎(chǔ)上采用非正交基的非正交化方式來對緊束縛（Tight Binding，TB）理論做修正，計(jì)算速度比傳統(tǒng)DFT 快兩三個(gè)數(shù)量級，核心是對哈密頓矩陣H的求解和對電子波函數(shù)λ的計(jì)算，能量表達(dá)式如下所示：

式中，occ 表示電子占據(jù)數(shù)，i表示具體的價(jià)電子表示原子核A 和B 之間的排斥能，kJ。電子波函數(shù)主要通過基函數(shù)的線性組合來求解：

式中，l是電子的空間坐標(biāo)，a=x，y，z表示空間維度。DFTB 方法采用經(jīng)驗(yàn)數(shù)值代替原本復(fù)雜的積分方程，大大降低了哈密頓矩陣的求解難度，可用于計(jì)算晶體結(jié)構(gòu)和密度［50］，但該方法沒有針對硼元素的經(jīng)驗(yàn)參數(shù)，無法模擬極端高溫高壓條件下的物理化學(xué)變化。J.R.Holden［51］用半經(jīng)驗(yàn)AM1（Austin Model Number 1）方法對堆積結(jié)構(gòu)進(jìn)行優(yōu)化，預(yù)測晶胞參數(shù)。

半經(jīng)驗(yàn)方法雖然可將計(jì)算規(guī)模增大，但還不能有效以及長時(shí)間模擬原子數(shù)量非常龐大的系統(tǒng)或者結(jié)構(gòu)更加靈活的體系，比如復(fù)雜分子表面勢能以及弱鍵能等；計(jì)算的可靠性與參數(shù)獲得的途徑密切相關(guān)，參數(shù)擬合大多依賴于DFT 計(jì)算；由于采用了一系列近似，應(yīng)用于含有強(qiáng)吸電子基團(tuán)（如硝基）的體系時(shí)，由于電荷密度不對稱，計(jì)算誤差會大大增加。

3.2 勢能面搜索預(yù)測晶體結(jié)構(gòu)

晶體結(jié)構(gòu)預(yù)測是在固定化學(xué)組分情況下確定最有可能的原子排列方式和晶胞參數(shù)，這一過程中用到自然界中占據(jù)統(tǒng)治地位的能量最低原理，即在固定組分的狀態(tài)子空間中搜索處于能量曲面極小值的排列方式。從理論上預(yù)測晶體結(jié)構(gòu)一直是化學(xué)、物理和材料科學(xué)中的一大難題，隨著理論方法和計(jì)算技術(shù)的完善，發(fā)展了許多晶體結(jié)構(gòu)預(yù)測的方法，為研究含能晶體的結(jié)構(gòu)及其性質(zhì)注入活力。2016 年，多個(gè)單位參加了對有機(jī)晶體結(jié)構(gòu)的第六次盲測［52］，分子結(jié)構(gòu)均來自第一原理計(jì)算和劍橋結(jié)構(gòu)數(shù)據(jù)庫，晶體結(jié)構(gòu)通過隨機(jī)搜索法、遺傳算法或模擬退火法產(chǎn)生。

3.2.1 模擬退火法

模擬退火法是S.Kirkpatrick 等［53］于1983 年基于金屬退火機(jī)理而建立的全局最優(yōu)化方法，通過對給定初始狀態(tài)進(jìn)行模擬加溫，生成一系列關(guān)聯(lián)的新解并評估與初始解的能量差，再根據(jù)Metropolis 接受準(zhǔn)則進(jìn)行結(jié)構(gòu)篩選［54］，得到一系列給定溫度下允許存在新解。例如在溫度T時(shí)，由當(dāng)前狀態(tài)i產(chǎn)生新狀態(tài)j，能量分別為Ei和Ej，若Ei＞Ej則接受新狀態(tài)為當(dāng)前狀態(tài)，否則，以概率p來接受狀態(tài)，其中k為Boltzmann 常數(shù)。

在模擬退火的開始階段通常以升高溫度來實(shí)現(xiàn)初始解在能量曲面上的上坡游走，克服勢壘走向周圍的局域最小值；這些走離初始解位置的新解通常都處于能量面中的坡上，需要進(jìn)一步通過降溫處理讓它們再次在能量曲面上進(jìn)行下陂游走，走向周圍的局域態(tài)。升溫和降溫時(shí)，溫度變化遵循如下關(guān)系式：

式中，Tnew和Told分別是變溫前后的溫度，K；Th和Tc分別是自定義的加熱因子和冷卻因子。

J.Yin［55］和G.Z.Zhao［56］分別用模擬退火技術(shù)預(yù)測了1′-氨基-1′H-1，5′-雙四唑（ATT）及其溴氰衍生物和2，4，6-三硝基-1，3，5-三氮-1，3，5-三氧化六元氮雜環(huán)（TNTATO）的晶體結(jié)構(gòu)。劉英哲等［57］基于B3LYP方法和6-31G（d）基組建立了適用于全氮結(jié)構(gòu)的力場，用Monte-Carlo 方法進(jìn)行空間堆積，預(yù)測晶體結(jié)構(gòu)并計(jì)算N4、N6、N8、N10、N12五種籠型全氮分子的晶體密度分別是1.81，2.08，2.47，2.46，2.57 g·cm-3。 J.F.Moxnes［58］用Polymorph Predictor 方法和COMPASS力場計(jì)算晶體密度并與實(shí)驗(yàn)數(shù)據(jù)比較（見表1），對大多數(shù)空間群，模擬退火法高估了密度，平均絕對誤差0.07～0.08 g·cm-3。

模擬退火法是早期含能晶體結(jié)構(gòu)預(yù)測常用方法，這種基于初始結(jié)構(gòu)演化的躍遷勢壘方法對人為選擇的初始結(jié)構(gòu)具有依賴性，考慮空間群不夠?qū)挿?，難以克服退火過程的高能量勢壘，需要進(jìn)行多次迭代直到所得結(jié)構(gòu)滿足所設(shè)置的收斂標(biāo)準(zhǔn)，在計(jì)算效率上有很大缺陷。

表1 Polymorph Predictor 方法計(jì)算的密度與實(shí)驗(yàn)數(shù)據(jù)對比Table 1 Comparison of the calculated densities by Polymorph Predictor methods and experimental data

3.2.2 隨機(jī)搜索法

隨機(jī)搜索方法通過在狀態(tài)空間盡可能多地進(jìn)行無區(qū)別隨機(jī)采樣，產(chǎn)生大量初始結(jié)構(gòu)Xst，參照（11）式在其周圍以一定的步長隨機(jī)搜索可行解：

式中，t為當(dāng)前迭代次數(shù)，rand( )為隨機(jī)數(shù)函數(shù)，a和b為隨機(jī)數(shù)的取值范圍；然后進(jìn)行精細(xì)結(jié)構(gòu)優(yōu)化使它們走向鄰近的能量極小值，經(jīng)過不斷迭代從而完成晶體結(jié)構(gòu)的預(yù)測［61］。

第一原理隨機(jī)搜索方法將DFT-D 和AIRSS結(jié)合，在晶體結(jié)構(gòu)預(yù)測中顯現(xiàn)出了強(qiáng)大的能力。如M.Zilka［62］用該方法成功的預(yù)測了有機(jī)分子m-氨基-苯甲酸（m-ABA）的晶體結(jié)構(gòu)，與X-ray 衍射和固態(tài)核磁共振的低溫實(shí)驗(yàn)結(jié)果吻合，并能區(qū)分不同類型氫鍵的作用。T.J.Lin［63］用AIRSS 方法發(fā)現(xiàn)了甲醇新的δ相候選結(jié)構(gòu)，通過DFT-D 計(jì)算確定了相變壓力和CH…O 氫鍵在高壓下晶體結(jié)構(gòu)的穩(wěn)定性中起主導(dǎo)作用。

隨機(jī)搜索方法的優(yōu)點(diǎn)是對目標(biāo)函數(shù)無特殊要求，通用性強(qiáng)。但計(jì)算需要足夠大的樣本方可獲得較精確的結(jié)構(gòu)。依據(jù)晶格能的計(jì)算得到晶體結(jié)構(gòu)之間的能量差很小，隨機(jī)搜索方法僅按照能量進(jìn)行篩選很容易漏掉其他熱力學(xué)穩(wěn)定結(jié)構(gòu)。在搜索進(jìn)程中增加自適應(yīng)功能，使算法能夠根據(jù)當(dāng)前解自動計(jì)算下次迭代搜索的方向和步長，可有效改進(jìn)搜索能力和精度。

3.2.3 進(jìn)化算法

進(jìn)化算法是通過不斷地產(chǎn)生候選晶體結(jié)構(gòu)種群，并經(jīng)過生物進(jìn)化理論中的遺傳繁衍、基因變異、基因重組和自然選擇等基本法則不斷更新和優(yōu)化這一種群，直至搜索到滿足要求的晶體結(jié)構(gòu)。在這一過程中，需要設(shè)定一個(gè)更新法則來決定如何進(jìn)行繁衍、變異和重組，同時(shí)還需要一個(gè)適合的選擇法則來決定個(gè)體的存活幾率、變異幾率和雜交幾率，以保證個(gè)體進(jìn)化能夠朝著優(yōu)化的方向進(jìn)行，該思想在A.R.Oganov 等開發(fā)的晶體結(jié)構(gòu)預(yù)測軟件USPEX（Universal Structure Predictor：Evolutionary Xtallography）中實(shí)現(xiàn)［64］。USPEX通過指紋函數(shù)描述晶體結(jié)構(gòu)：

式中，Z為原子相對質(zhì)量；Rij為兩原子間距離，nm；V為單位晶胞體積，nm3；N 為單位晶胞內(nèi)所含有的原子數(shù)，R是一個(gè)可變參數(shù)。不同結(jié)構(gòu)之間的相似性可由指紋函數(shù)空間定義的距離判斷：

USPEX 軟件基于第一性原理和進(jìn)化算法，僅需提供材料的化學(xué)組分就可以得到能量、體積、硬度、介電常數(shù)、能帶寬度和磁矩等性質(zhì)，已經(jīng)預(yù)測了很多無機(jī)晶體結(jié)構(gòu)［65-68］，但用于含能晶體預(yù)測的報(bào)道較少。王洪波［69］用該方法預(yù)測了含氮高能量高密度材料的晶體結(jié)構(gòu)，預(yù)言了SiC2N4和Si2CN4的高致密相。

3.2.4 粒子群優(yōu)化算法

粒子群優(yōu)化算法是在進(jìn)化算法的基礎(chǔ)上發(fā)展起來的一種基于種群搜索策略的全局優(yōu)化算法，它和模擬退火算法相似，都以隨機(jī)解為起點(diǎn)，通過迭代手段搜索最優(yōu)解，在該過程中需要定義適應(yīng)度來標(biāo)記和評估嘗試解的品質(zhì)。PSO 不需要進(jìn)行“交叉”和“變異”操作，而直接通過當(dāng)前搜索的最值來確定下一步的搜索方向，達(dá)到對全局最優(yōu)解的搜索［70］。該算法中粒子更新的速度和位置由以下公式?jīng)Q定：

在結(jié)構(gòu)搜索領(lǐng)域，一些全局性參數(shù)優(yōu)化方法與考慮對稱性的結(jié)構(gòu)產(chǎn)生技術(shù)、判斷相似結(jié)構(gòu)的幾何結(jié)構(gòu)因子的引入、以及提高群體多樣性的技術(shù)等相結(jié)合，有效的克服了前期晶體結(jié)構(gòu)預(yù)測中難以保持種群多樣性的問題，可以保證對勢能面進(jìn)行有效的探索。后來發(fā)展的基于高斯過程回歸的機(jī)器學(xué)習(xí)勢函數(shù)，將群體智能與機(jī)器學(xué)習(xí)方法的結(jié)合，能進(jìn)行大尺寸復(fù)雜材料的結(jié)構(gòu)預(yù)測［73］。

4 基于數(shù)學(xué)統(tǒng)計(jì)方法預(yù)測密度

4.1 定量構(gòu)效關(guān)系模型

通過數(shù)學(xué)和統(tǒng)計(jì)學(xué)手段找出含能材料的化學(xué)結(jié)構(gòu)與密度等性質(zhì)之間的量變規(guī)律，或得到構(gòu)效關(guān)系的數(shù)學(xué)方程，可以為密度預(yù)測提供理論依據(jù)［74-76］，核心問題是采用有效算法建立構(gòu)效關(guān)系模型。本文介紹已經(jīng)用于含能材料密度預(yù)測模型的建立方法：回歸分析、人工神經(jīng)網(wǎng)絡(luò)和改進(jìn)的遺傳算法。

4.1.1 回歸分析

回歸分析是處理變量之間相關(guān)關(guān)系的常用數(shù)理統(tǒng)計(jì)方法［77］，可通過后者的已知或設(shè)定值，去估計(jì)或預(yù)測前者的均值。在回歸分析中，如果有兩個(gè)或兩個(gè)以上的自變量，就稱為多元回歸［78］。來蔚鵬等［79］對30 種芳香系單質(zhì)炸藥的偶極矩、最高占據(jù)軌道能量、分子總能量等8 個(gè)描述符進(jìn)行多元線性回歸計(jì)算，建立QSPR 模型預(yù)測密度，訓(xùn)練集和測試集的預(yù)測值與實(shí)驗(yàn)值的平均誤差分別為3.33%和2.94%。

支持向量機(jī)（Support Vector Machine，SVM）是根據(jù)有限的樣本信息，在模型復(fù)雜性和學(xué)習(xí)能力之間尋求最佳折中，以獲得良好的推廣能力。支持向量回歸（Support Vector Regression，SVR）是用于解決回歸問題的支持向量機(jī)。宗朝霞等［74］基于遺傳算法的變量篩選和支持向量機(jī)預(yù)測呋咱系含能化合物的密度，隨機(jī)選取85%的待測呋咱系含能化合物作為訓(xùn)練集，其余為測試集，預(yù)測結(jié)果平均相對誤差分別為1.16%和2.12%。S.N.Hwang［80］基于多元線性回歸和支持向量機(jī)建立QSPR 模型預(yù)測了3609 種含能材料的密度，結(jié)果與Ammo 的基團(tuán)貢獻(xiàn)法相當(dāng)。

支持向量機(jī)和支持向量回歸是目前機(jī)器學(xué)習(xí)領(lǐng)域用得較多的方法，序列最小優(yōu)化算法（Sequential Minimal Optimization，SMO）的提出有效的解決了支持向量機(jī)方法實(shí)現(xiàn)復(fù)雜、效率低等問題，特別在處理大數(shù)據(jù)量的實(shí)際問題時(shí)體現(xiàn)出了較好的精度和運(yùn)算效率，可用于含能材料設(shè)計(jì)。

4.1.2 人工神經(jīng)網(wǎng)絡(luò)

人工神經(jīng)網(wǎng)絡(luò)是人工智能研究一種數(shù)學(xué)模型［81］，在解決復(fù)雜環(huán)境系統(tǒng)非線性關(guān)系模擬等眾多問題上展示出強(qiáng)大的功能，越來越多的應(yīng)用于炸藥性能的研究［82-84］。K.Tatyana 等［85］運(yùn)用人工智能和模式識別方法，通過結(jié)構(gòu)基團(tuán)分析與識別，從大量已知結(jié)構(gòu)數(shù)據(jù)庫得到規(guī)律，對含能材料的性能進(jìn)行初步預(yù)測。人工神經(jīng)網(wǎng)絡(luò)可以正確反映配方組成和溫度對裝藥密度的影響規(guī)律，對裝藥密度進(jìn)行預(yù)測［86］。王國棟等［87-88］結(jié)合DFT 與人工神經(jīng)網(wǎng)絡(luò)，以平均極化率、平均四極矩和偶極矩為描述符，構(gòu)建炸藥密度的預(yù)測模型，預(yù)測16 種炸藥的密度與文獻(xiàn)值相對誤差為0.47%～6.6%。

神經(jīng)網(wǎng)絡(luò)是一個(gè)可以自學(xué)習(xí)并能夠“記憶”樣本所含信息的理論方法，得到廣泛的應(yīng)用，但它也存在自身的限制與不足，如易形成局部極小值而得不到全局最優(yōu)解［89］；訓(xùn)練次數(shù)多使得學(xué)習(xí)效率降低，收斂速率慢；隱節(jié)點(diǎn)的選取缺乏理論指導(dǎo)等。針對上述問題，國內(nèi)外已提出不少有效的改進(jìn)方法，較常用的有增加動量項(xiàng)、自適應(yīng)調(diào)節(jié)學(xué)習(xí)速率、引入陡度因子等。

4.1.3 改進(jìn)的遺傳算法

遺傳算法具有自組織、自適應(yīng)和學(xué)習(xí)性等特點(diǎn)，搜索過程以目標(biāo)函數(shù)值為評價(jià)標(biāo)準(zhǔn)，不受優(yōu)化函數(shù)連續(xù)性的約束。針對遺傳算法收斂速度慢的缺點(diǎn)，程娥等［90］提出一種偽并行、貪婪的、最優(yōu)解保存與自適應(yīng)參數(shù)調(diào)整相結(jié)合的改進(jìn)遺傳算法，對標(biāo)準(zhǔn)化后的399 個(gè)炸藥分子結(jié)構(gòu)描述符和爆轟數(shù)據(jù)進(jìn)行變量的選擇并建立定量“分子結(jié)構(gòu)-爆炸性能”關(guān)系（QSDR）模型（式（16）），從結(jié)構(gòu)描述符中篩選出13 個(gè)與密度相關(guān)性最大的參數(shù)，呋咱類化合物密度的預(yù)測值相對誤差在±2%以內(nèi)。

偽并行是將群體分成N 個(gè)子群，各自按照不同的交叉率和變異率進(jìn)行優(yōu)化，最后比較各個(gè)子群的最優(yōu)解得到整個(gè)群體的最優(yōu)解；最優(yōu)解保存策略是每次從所有子群中找出一個(gè)最優(yōu)個(gè)體，再將此個(gè)體注入每個(gè)子群替代各子群中最差個(gè)體，防止隨機(jī)選擇時(shí)適應(yīng)值最好個(gè)體未被選中的情況；自適應(yīng)參數(shù)調(diào)整即交叉率和變異率通過個(gè)體本身適應(yīng)度和種群整體性能的比較確定；貪婪策略通過對父代、子代4 個(gè)個(gè)體適應(yīng)度作比較，選擇適應(yīng)度較大的2 個(gè)作為子代個(gè)體，確保最佳個(gè)體被遺傳到下一代而不受破壞。改進(jìn)的遺傳算法與傳統(tǒng)遺傳算法相比，收斂速度快、尋優(yōu)能力強(qiáng)，誤差率小，且性能穩(wěn)定，有望在含能化合物研發(fā)中發(fā)揮更大作用。

4.2 經(jīng)驗(yàn)公式估算法

目前，國際上由經(jīng)驗(yàn)公式估算含能材料晶體密度的研究中，最有代表性和被廣泛引用的報(bào)道來自于M.H.Keshavarz。2007 年Keshavarz［91-92］由元素組成預(yù)測部分含能化合物的晶體密度，所用公式如下：

MW為分子摩爾質(zhì)量，g·mol-1；a、b、c、d、ni是C、H、N、O 原子及官能團(tuán)的數(shù)量；經(jīng)過擬合得到：

對不同的分子結(jié)構(gòu)、不同的官能團(tuán)及連接方式，還需要分別對上式進(jìn)行復(fù)雜的修正才能用于密度預(yù)測，例如對于分子結(jié)構(gòu)中一個(gè)C 原子連接三個(gè)硝基的情況，密度計(jì)算公式如下：

其余情況用ρ0，corr= 0.1233 + 0.8373ρ0修正。

2009 年，Keshavarz［93］對前期建立的密度計(jì)算公式修正如下：

用多元線性回歸方法從實(shí)驗(yàn)數(shù)據(jù)中尋找可調(diào)參數(shù)z1～z6，EI和ED是與特殊結(jié)構(gòu)相對應(yīng)的參數(shù)，可增加或降低晶體密度，得到密度計(jì)算公式：

將72 種化合物的密度預(yù)測結(jié)果與前期的經(jīng)驗(yàn)公式法和基團(tuán)加和法進(jìn)行比較，準(zhǔn)確度有較大提升。

2013 年，Keshavarz［94-97］在對元素組成與特定分子進(jìn)行修正的同時(shí)，綜合了定量構(gòu)效關(guān)系方法，考慮分子間的相互作用，應(yīng)用多元線性回歸模型，對含能化合物的密度進(jìn)行預(yù)測，方法的適用性和準(zhǔn)確度有較大提高。D.Frem［98］用式（23）的經(jīng)驗(yàn)關(guān)系方法預(yù)測的晶體密度值偏高。

CPG和CNG分別是特殊官能團(tuán)結(jié)構(gòu)對密度的正、負(fù)貢獻(xiàn)。

2015 年，Keshavarz［99］改進(jìn)了上述經(jīng)驗(yàn)關(guān)系法，將密度計(jì)算公式修正為：

式中，IMP 和DMP 分別是分子間相互作用對密度的修正，體現(xiàn)了特定分子中H 原子數(shù)nH和N 原子數(shù)nN的偏差，對不同的結(jié)構(gòu)取不同參數(shù)，參數(shù)的確定來源于實(shí)驗(yàn)數(shù)據(jù)。對172 種不同種類的含能化合物計(jì)算得到的均方根誤差為2.16%～11.0%。

2017 年，Keshavarz 發(fā)布了EMDB（Energetic Materials Designing Bench）軟件［100］，可以計(jì)算含能材料的30 多種物理化學(xué)性質(zhì)及爆轟性能。

2018 年，V.D.Ghule［101］課題組考慮了含能分子中含有兩個(gè)及以上NH2或NH3+基團(tuán)對密度的貢獻(xiàn)，將密度計(jì)算公式修正為如下形式：

式中，VC、VH、VN、VO是C、H、N、O 的原子體積，nm3。對119 個(gè)含苯環(huán)的含能材料密度進(jìn)行計(jì)算并與實(shí)驗(yàn)數(shù)據(jù)對比，其中27 個(gè)絕對誤差小于0.05 g·cm-3，硝基苯及其衍生物中54 個(gè)誤差范圍在3%以內(nèi)，12 個(gè)誤差范圍是3%～4%，8 個(gè)誤差范圍是5%～8%，離子鹽和共晶中25 個(gè)誤差在3%以內(nèi)，13 個(gè)誤差范圍是3%～4%，7個(gè)誤差在5%～7.4%。

QSPR 預(yù)測含能晶體密度時(shí)，訓(xùn)練樣本越多預(yù)測效果越好，而實(shí)際情況是很多炸藥的實(shí)驗(yàn)數(shù)據(jù)缺乏，影響了模型的精度［102-103］，很難明確給出方程的物理意義和分子間相互作用的模式及化學(xué)內(nèi)涵；經(jīng)驗(yàn)公式方法雖有直接、快速的特點(diǎn)，但普適性差，且忽略了分子間相互作用對晶體堆積結(jié)構(gòu)的影響，在原理上不夠準(zhǔn)確。

5 展望

綜上所述，每種密度計(jì)算方法在擁有諸多優(yōu)勢的同時(shí)都存在一定的局限性，如基于分子體積的方法在CHNO 類含能材料中應(yīng)用比較廣泛，密度計(jì)算的準(zhǔn)確度取決于能否正確的描述分子間和分子內(nèi)相互作用、空間位阻效應(yīng)、復(fù)雜的環(huán)結(jié)構(gòu)以及致爆基，如NO2、NH2、OH、N3、ONO2、NHNO2等對密度的影響。MSEP 優(yōu)于IED 的計(jì)算結(jié)果，一般會略微低估密度，基組選擇對密度的計(jì)算結(jié)果影響較小?；诹龅娜S晶體堆積方法預(yù)測的密度一般會比實(shí)驗(yàn)值略高，晶體結(jié)構(gòu)搜索要面對巨大的狀態(tài)空間和高度復(fù)雜的能量曲面，發(fā)展高效可靠的晶體結(jié)構(gòu)預(yù)測方法一直是亟待解決的問題。經(jīng)驗(yàn)公式法雖然直接但對新發(fā)現(xiàn)的分子或某些特殊結(jié)構(gòu)的適用性需驗(yàn)證，QSPR 預(yù)測含能晶體密度的模型精度需要提高。以下幾種方案結(jié)合量子物理、計(jì)算化學(xué)、人工智能三者的優(yōu)勢，可顯著提高含能材料密度計(jì)算的效率和可靠性，實(shí)現(xiàn)含能材料的高通量設(shè)計(jì)：

（1）將群體智能搜索與密度泛函理論的色散修正和機(jī)器學(xué)習(xí)勢函數(shù)方法結(jié)合，進(jìn)行含能晶體的結(jié)構(gòu)和密度預(yù)測：先通過高斯過程回歸得到勢函數(shù)，用于勢能面快速搜索，提高效率；采用全空間密網(wǎng)格初篩保證初篩樣本的完備性，采用剛性分子堆積技術(shù)或晶體配位搜索技術(shù)減小自由度，并適用于多種成鍵類型，提高成功率；通過自由能表面全局最小化搜索得到備選晶體結(jié)構(gòu)后，再用DFT-D 方法對晶體結(jié)構(gòu)和密度精細(xì)計(jì)算，提高置信度。

（2）定量構(gòu)效關(guān)系與人工智能方法結(jié)合，可以方便準(zhǔn)確的預(yù)測含能晶體密度：建立神經(jīng)網(wǎng)絡(luò)模型描述復(fù)雜的相關(guān)關(guān)系，結(jié)合數(shù)據(jù)挖掘，可以進(jìn)一步提高模型的代表性和預(yù)測精度；該方案需要高質(zhì)量的含能材料結(jié)構(gòu)與性能的實(shí)驗(yàn)數(shù)據(jù)做支持，為分子設(shè)計(jì)和結(jié)構(gòu)-性能關(guān)系的預(yù)測提供訓(xùn)練依據(jù)。

（3）第一原理和機(jī)器學(xué)習(xí)結(jié)合（QMML 方法）：通過高通量集成計(jì)算產(chǎn)生大批量的含能分子及復(fù)合體系，將相關(guān)結(jié)構(gòu)和性質(zhì)的計(jì)算結(jié)果形成數(shù)據(jù)庫，用于分子表面靜電勢修正，基于機(jī)器學(xué)習(xí)動態(tài)更新樣本結(jié)構(gòu)，一方面通過開發(fā)專用基函數(shù)和高精度贗勢，提高第一原理大規(guī)模并行計(jì)算能力，進(jìn)行結(jié)構(gòu)構(gòu)造和性能預(yù)估；另一方面，高通量計(jì)算與數(shù)據(jù)挖掘融合，可能發(fā)現(xiàn)隱藏的材料“基因”-從原子排列到相結(jié)構(gòu)再到顯微組織形成，最后到材料宏觀性能與使用壽命之間的關(guān)系。