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

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

1 引言

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

2 基于分子體積計算密度

2.1 等電子密度面法

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

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

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

2.2 分子表面靜電勢校正法

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

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

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

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

2.3 基團加和法

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

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

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

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

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

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

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

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

經(jīng)過DFT-D 校正后，靜電主導(dǎo)的弱相互作用，如氫鍵、鹵鍵等問題，計算精度有較大改進，已用于計算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ù)計算結(jié)果與實驗的誤差通常在2%以內(nèi)［31］。L.Zhang等［32］采用PBE交換關(guān)聯(lián)勢計算HMX 的α、β、δ和γ相的晶體結(jié)構(gòu)，體積誤差分別為-2.36%，-1.00%，+0.12%和-0.3%，采用vdW-D3方法計算了53 種炸藥晶體結(jié)構(gòu)，體積誤差＜±4%［33］。H.Lin［34］、H.C.S.Chan［35］、P.Panini［36］也采用DFT-D對含能晶體或共晶進行結(jié)構(gòu)和密度預(yù)測。

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

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

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

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

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

3.1.3 半經(jīng)驗方法

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

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

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

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

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

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

3.2.1 模擬退火法

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

在模擬退火的開始階段通常以升高溫度來實現(xiàn)初始解在能量曲面上的上坡游走，克服勢壘走向周圍的局域最小值；這些走離初始解位置的新解通常都處于能量面中的坡上，需要進一步通過降溫處理讓它們再次在能量曲面上進行下陂游走，走向周圍的局域態(tài)。升溫和降溫時，溫度變化遵循如下關(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 方法進行空間堆積，預(yù)測晶體結(jié)構(gòu)并計算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力場計算晶體密度并與實驗數(shù)據(jù)比較（見表1），對大多數(shù)空間群，模擬退火法高估了密度，平均絕對誤差0.07～0.08 g·cm-3。

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

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

3.2.2 隨機搜索法

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

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

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

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

3.2.3 進化算法

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

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

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

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

粒子群優(yōu)化算法是在進化算法的基礎(chǔ)上發(fā)展起來的一種基于種群搜索策略的全局優(yōu)化算法，它和模擬退火算法相似，都以隨機解為起點，通過迭代手段搜索最優(yōu)解，在該過程中需要定義適應(yīng)度來標記和評估嘗試解的品質(zhì)。PSO 不需要進行“交叉”和“變異”操作，而直接通過當前搜索的最值來確定下一步的搜索方向，達到對全局最優(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ù)測中難以保持種群多樣性的問題，可以保證對勢能面進行有效的探索。后來發(fā)展的基于高斯過程回歸的機器學(xué)習勢函數(shù)，將群體智能與機器學(xué)習方法的結(jié)合，能進行大尺寸復(fù)雜材料的結(jié)構(gòu)預(yù)測［73］。

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

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

通過數(shù)學(xué)和統(tǒng)計學(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ò)和改進的遺傳算法。

4.1.1 回歸分析

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

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

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

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

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

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

4.1.3 改進的遺傳算法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

5 展望

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

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

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

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