含能晶体密度预测的研究进展

王丽莉，熊 鹰，谢炜宇，牛亮亮，张朝阳

（1.中国工程物理研究院计算机应用研究所，四川 绵阳 621999；2.中国工程物理研究院化工材料研究所，四川 绵阳 621999）

1 引言

含能材料密度的理论预测方法可分成以下三大类：一是基于分子体积计算密度，由于分子体积没有明确的物理定义，因此存在多种计算方法，如等电子密度面（Isosurface of Electron Density，IED）方法通常将分子摩尔体积定义为0.00l a.u.电子密度等值面轮廓范围内的体积［1］，分子表面静电势（Molecular Surface Electrostatic Potentials，MESP）校正法是在密度求解公式中加入通过实验数据拟合的系数表征分子间相互作用，基团加和（Group Additivity Methods，GAM）法是由原子和官能团体积的贡献估算分子体积等，其中计算电子密度的方法多数采用密度泛函理论（Density Functional Theory，DFT）、二阶微扰（Second-order Moller-Plesset Perturbation，MP2）理论、耦合簇（Coupled Cluster with Singles and Doubles，CCSD）理论等。二是基于晶体体积计算密度，晶体结构预测（Crystal Structure Prediction，CSP）结合了能量和结构的精确计算方法与全局结构搜索方法，其中结构与能量的计算方法包括密度泛函理论的色散校正（Dispersion Correction Density Functional Theory，DFT-D），分子力场与分子动力学方法（Molecular Dynamics，MD），半经验分子轨道方法（Semi-empirical Molecular Orbital，MO）等；结构搜索的代表性算法有模拟退火法（Simulated Annealing，SA）、随机结构搜索法（Ab Initio Random Structure Searching，AIRSS）、进化算法（Evolutionary Approach，EA）和粒子群优化算法（Particle Swarm Optimization，PSO）等。三是基于数学和统计方法估算含能材料密度，如通过多元线性回归（Multiple Linear Regressin，MLR）、人工神经网络（Artificial Neural Networks，ANN）和改进的遗传算法（Genetic Algorithm，GA）等建立快速稳定的定量构效关系（Quantitative Structure Property Relationship，QSPR）模型，或根据经验公式由元素组成和官能团数目估算晶体密度等。随着云计算的普及和大数据时代的到来，机器学习（Machine Learning，ML）及人工智能（Artificial Intelligence，AI）技术将推动含能晶体结构、性质以及感度之间关系的研究，并为开发新型钝感高能炸药提供依据。

2 基于分子体积计算密度

2.1 等电子密度面法

等电子密度面法计算密度的主要步骤如下：优化含能分子确保其处于能量最低点，然后计算振动光谱和相应热力学信息评估结构稳定性［2］，用Monte-Carlo方法计算每个分子的平均摩尔体积V（0.001），理论密度用分子摩尔质量M与V（0.001）之比求得［3］。该方法对分子形状和体积的描述依赖电子密度值的选择，对于处于凝聚态的分子，考虑到分子间压积会使范德华表面有一定穿透，也有用0.002 a.u.或0.0015 a.u.的电子密度等值面作为范德华表面［4-6］。

密度泛函理论的B3LYP 方法和6-31G**基组是等电子密度面法计算分子晶体密度的主要方法［7］，准确性优于半经验MO 方法［8］，经济性优于MP2 和CCSD（T）等高水平电子密度方法［9］。对CHNO 型分子晶体，该方法计算中性分子的理论密度与实验值的平均误差和均方根误差分别为1%和3.7%，分子中含氮量超过50%或含氮环的误差为-1.8%和3.4%，离子晶体的误差为-5.1%和6.6%，通过对体积公式进行修正，离子晶体平均误差和均方根误差减小至-1.1%和4.8%［10］。

等电子密度面法没有考虑分子间空隙及各种复杂的相互作用，用于密度这种宏观性质计算，原理受质疑。由于该方法对极性分子的计算误差偏大，后期根据经验进行参数的计算拟合，对不同的分子结构进行体积修正，精确度有所提高。

2.2 分子表面静电势校正法

P.Politzer 等［11-12］提出修正分子表面静电势参数法，在密度求解公式中分别考虑了与分子表面正、负电荷分布情况相关的分子表面静电势平衡系数v和总方差计算公式见式（1）：

式中，M为分子摩尔质量，g·mol-1；α、β、γ是衰减系数，由计算拟合得到。该方法预测密度与实验值误差一般小于0.05 g·cm-3，被广泛应用于CHNO 型化合物密度的预测［13］。如H.Singh［14］、P.Ma［15］和Y.J.Shu等［16］用Politzer 方法计算了2，4，6-三硝基甲苯（TNT）、环三亚甲基三硝胺（RDX）、环四亚甲基四硝胺（HMX）等多种含能分子的结构和密度，都与实验数据吻合较好。2013 年，B.M.Rice［17］用Politzer 方法对前期工作进行改进［10］，并将计算结果分组与实验数据比较，其中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 修正更接近实验值，多数误差小于3%，硝酸脂和包含笼状或刚性环结构的误差较大。

分子表面静电势校正方法是目前计算含能材料密度的主要方法，计算精度优于等电子密度面法和基团加和法，其可靠性己经被证明。

2.3 基团加和法

基团加和法用组成分子的原子和官能团的标准体积之和估算分子体积，再用分子摩尔质量求密度。H.L.Ammon 等［21］从11555 个晶体结构数据中提取了78 种原子和官能团的体积，建立标准体积数据库，用于计算CHNOF 类化合物的晶体密度，后来又扩充和升级该数据库［22-23］，用于更多有机物的密度预测。D.Mathieu［24］以Ammon 数据库为基础，用基团加和法预测密度与实验数据的平均误差接近2%，只有部分密堆积结构的估算密度误差达到6%～8%。S.Beaucamp［25-26］用带电基团体积加和法估算离子盐和水合物的晶体密度，平均误差小于2.5%，通过修正氢键和环状结构等对分子体积的贡献［27］，42880 个晶体估算密度的平均误差为2%。

基团加和法中体积的贡献值是从大量实验数据中总结得出，缺乏实验数据难以获得准确的经验参数。不仅如此，该方法不考虑分子构型和分子间相互作用，无法区分同分异构体的密度差异，忽略晶型及温度对密度的影响，对非电中性含能化合物，含硅、硼等元素的化合物，或部分密堆积结构含能晶体的计算偏差较大，估算范围有一定限制。

3 基于晶体结构计算密度

3.1 量化方法优化晶体结构

从理论上预测含能晶体的结构，可以得到密度。晶体结构的求解通常需要高精度量化方法进行结构优化来实现局域最小值精修和能量稳定性计算。

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

密度泛函理论对分子间弱相互作用的描述不完全准确，为了能够精确的计算含能分子晶体的物性，需要加入分子间的色散校正。非局域性色散校正、London 色散校正、球原子模型等是对DFT 的有效修正。其中DFT-D 方法在量化程序中得到广泛支持，通常用零阻尼函数（式（5））来调节色散校正在近程、中程距离时的行为，避免近距离时校正能（式（4））较大导致的重复问题［28］。

其中RAB代表AB 原子间的距离，sn是截断半径，nm；C是原子间色散校正系数；sn和sr，n是刻度因子；γ是预设常数。随原子间距离减小，f逐渐为0，故称零阻尼。

经过DFT-D 校正后，静电主导的弱相互作用，如氢键、卤键等问题，计算精度有较大改进，已用于计算5，5′-联四唑-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］，晶胞参数计算结果与实验的误差通常在2%以内［31］。L.Zhang等［32］采用PBE交换关联势计算HMX 的α、β、δ和γ相的晶体结构，体积误差分别为-2.36%，-1.00%，+0.12%和-0.3%，采用vdW-D3方法计算了53 种炸药晶体结构，体积误差＜±4%［33］。H.Lin［34］、H.C.S.Chan［35］、P.Panini［36］也采用DFT-D对含能晶体或共晶进行结构和密度预测。

DFT-D 是计算含能材料晶体密度最适合的DFT方法，计算精度由泛函和色散校正模型的优良性、两者对色散作用描述的互补性两个因素决定，误差明显小于传统的DFT 方法。DFT-D 的可行性与精度经过了系统和全面的验证，但计算量随体系的增大呈三次方或二次方增长，比较耗时。采用对称适应微扰理论（Symmetry Adapted Perturbation Theory，SAPT）把相互作用能直接分解为静电相互作用、诱导能和色散能、对称性微扰理论定义的交换能，可以计算更大规模的原子模型。R.Podeszwa［37］采用基于DFT 单体描述的SAPT 准确预测了RDX 的晶体结构，计算了FOX-7 的完整势能面。

3.1.2 分子动力学结合第一原理

提高DFT 方法计算规模和效率的方法之一是与分子力场和分子动力学方法结合。MD 方法可研究单质炸药及其复合体系的结构、晶体生长机制［38］、能量、热稳定性、堆积密度、爆轰性能、溶剂对晶体的影响、力学性能及其温度效应，并与感度关联研究其安全性能，计算结果与实验数据越来越接近［39］。

分子动力学结合第一原理计算在晶体密度求解方面取得了极大的进展。P.Srinivasan［40］用B3LYP 方法和6-311G*基组优化2，4-二硝基苯甲酸（DNBA）的分子结构，再用MOLPAK（MOLecular PAcKing）建立空间堆积模型，进行UMD 力场晶格能最小化［41］，晶体结构预测结果与X-ray 衍射实验结果的偏差很小，密度的计算结果与实验值误差为2.98%。D.N.Kaushik［42］用第一原理处理有机晶体中的单分子和短程二体相互作用，用极化分子力场处理多体和长程相互作用，建立了基于片段的杂化多体相互作用模型，晶胞参数计算结果优于色散修正B3LYP-D*和AMOEBA 极化力场方法。J.S.Li［43］在B3LYP 方法和6-31++G**基组计算分子结构的基础上，用Dreiding 力场和模拟退火法预测了2，4，6，8-四硝基-1，3，5，7-四氮杂环丁烷（TNTAzC）和2，3，5，6，7，8-六硝基-1，4-二氮杂环丁烷（HNDAzC）的晶体密度分别为2.156 g·cm-3和1.975 g·cm-3。力场选择对晶体堆积和密度计算非常关键，Dreiding 力场对84%的不含芳香环结构的密度预测值与实验值的偏差在5%以内，PCFF 力场对83%的含芳香环结构的密度预测值与实验值的偏差在4%以内［44］，广义Amber 力场计算FOX-7、RDX、CL-20、HMX 等的密度比实验值低5%～8%［45］，CVFF 力场和DMACRYS 力场也用于含芳香环结构的含能晶体预测，密度误差小于7%［46-47］。

分子动力学的核心在于求解正则运动方程，分子内和分子间的作用力都包含在分子力场中，因此计算的准确性依赖所使用的力场［48-49］。如果缺乏完全适用的力场参数，无法准确描述分子间的相互作用，致使晶体结构预测存在较大偏差，要提高晶体密度计算的准确性需自研发经验势。

3.1.3 半经验方法

半经验方法通过忽略分子积分计算或引进经验参量对Hartree-Fock 方程的近似求解，如电荷自恰密度泛函紧束缚（Self-consistent Charge Density Functional Tight-binding，SCC-DFTB）方法是在DFT 基础上采用非正交基的非正交化方式来对紧束缚（Tight Binding，TB）理论做修正，计算速度比传统DFT 快两三个数量级，核心是对哈密顿矩阵H的求解和对电子波函数λ的计算，能量表达式如下所示：

式中，occ 表示电子占据数，i表示具体的价电子表示原子核A 和B 之间的排斥能，kJ。电子波函数主要通过基函数的线性组合来求解：

式中，l是电子的空间坐标，a=x，y，z表示空间维度。DFTB 方法采用经验数值代替原本复杂的积分方程，大大降低了哈密顿矩阵的求解难度，可用于计算晶体结构和密度［50］，但该方法没有针对硼元素的经验参数，无法模拟极端高温高压条件下的物理化学变化。J.R.Holden［51］用半经验AM1（Austin Model Number 1）方法对堆积结构进行优化，预测晶胞参数。

半经验方法虽然可将计算规模增大，但还不能有效以及长时间模拟原子数量非常庞大的系统或者结构更加灵活的体系，比如复杂分子表面势能以及弱键能等；计算的可靠性与参数获得的途径密切相关，参数拟合大多依赖于DFT 计算；由于采用了一系列近似，应用于含有强吸电子基团（如硝基）的体系时，由于电荷密度不对称，计算误差会大大增加。

3.2 势能面搜索预测晶体结构

晶体结构预测是在固定化学组分情况下确定最有可能的原子排列方式和晶胞参数，这一过程中用到自然界中占据统治地位的能量最低原理，即在固定组分的状态子空间中搜索处于能量曲面极小值的排列方式。从理论上预测晶体结构一直是化学、物理和材料科学中的一大难题，随着理论方法和计算技术的完善，发展了许多晶体结构预测的方法，为研究含能晶体的结构及其性质注入活力。2016 年，多个单位参加了对有机晶体结构的第六次盲测［52］，分子结构均来自第一原理计算和剑桥结构数据库，晶体结构通过随机搜索法、遗传算法或模拟退火法产生。

3.2.1 模拟退火法

模拟退火法是S.Kirkpatrick 等［53］于1983 年基于金属退火机理而建立的全局最优化方法，通过对给定初始状态进行模拟加温，生成一系列关联的新解并评估与初始解的能量差，再根据Metropolis 接受准则进行结构筛选［54］，得到一系列给定温度下允许存在新解。例如在温度T时，由当前状态i产生新状态j，能量分别为Ei和Ej，若Ei＞Ej则接受新状态为当前状态，否则，以概率p来接受状态，其中k为Boltzmann 常数。

在模拟退火的开始阶段通常以升高温度来实现初始解在能量曲面上的上坡游走，克服势垒走向周围的局域最小值；这些走离初始解位置的新解通常都处于能量面中的坡上，需要进一步通过降温处理让它们再次在能量曲面上进行下陂游走，走向周围的局域态。升温和降温时，温度变化遵循如下关系式：

式中，Tnew和Told分别是变温前后的温度，K；Th和Tc分别是自定义的加热因子和冷却因子。

J.Yin［55］和G.Z.Zhao［56］分别用模拟退火技术预测了1′-氨基-1′H-1，5′-双四唑（ATT）及其溴氰衍生物和2，4，6-三硝基-1，3，5-三氮-1，3，5-三氧化六元氮杂环（TNTATO）的晶体结构。刘英哲等［57］基于B3LYP方法和6-31G（d）基组建立了适用于全氮结构的力场，用Monte-Carlo 方法进行空间堆积，预测晶体结构并计算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力场计算晶体密度并与实验数据比较（见表1），对大多数空间群，模拟退火法高估了密度，平均绝对误差0.07～0.08 g·cm-3。

模拟退火法是早期含能晶体结构预测常用方法，这种基于初始结构演化的跃迁势垒方法对人为选择的初始结构具有依赖性，考虑空间群不够宽泛，难以克服退火过程的高能量势垒，需要进行多次迭代直到所得结构满足所设置的收敛标准，在计算效率上有很大缺陷。

表1 Polymorph Predictor 方法计算的密度与实验数据对比Table 1 Comparison of the calculated densities by Polymorph Predictor methods and experimental data

3.2.2 随机搜索法

随机搜索方法通过在状态空间尽可能多地进行无区别随机采样，产生大量初始结构Xst，参照（11）式在其周围以一定的步长随机搜索可行解：

式中，t为当前迭代次数，rand( )为随机数函数，a和b为随机数的取值范围；然后进行精细结构优化使它们走向邻近的能量极小值，经过不断迭代从而完成晶体结构的预测［61］。

第一原理随机搜索方法将DFT-D 和AIRSS结合，在晶体结构预测中显现出了强大的能力。如M.Zilka［62］用该方法成功的预测了有机分子m-氨基-苯甲酸（m-ABA）的晶体结构，与X-ray 衍射和固态核磁共振的低温实验结果吻合，并能区分不同类型氢键的作用。T.J.Lin［63］用AIRSS 方法发现了甲醇新的δ相候选结构，通过DFT-D 计算确定了相变压力和CH…O 氢键在高压下晶体结构的稳定性中起主导作用。

随机搜索方法的优点是对目标函数无特殊要求，通用性强。但计算需要足够大的样本方可获得较精确的结构。依据晶格能的计算得到晶体结构之间的能量差很小，随机搜索方法仅按照能量进行筛选很容易漏掉其他热力学稳定结构。在搜索进程中增加自适应功能，使算法能够根据当前解自动计算下次迭代搜索的方向和步长，可有效改进搜索能力和精度。

3.2.3 进化算法

进化算法是通过不断地产生候选晶体结构种群，并经过生物进化理论中的遗传繁衍、基因变异、基因重组和自然选择等基本法则不断更新和优化这一种群，直至搜索到满足要求的晶体结构。在这一过程中，需要设定一个更新法则来决定如何进行繁衍、变异和重组，同时还需要一个适合的选择法则来决定个体的存活几率、变异几率和杂交几率，以保证个体进化能够朝着优化的方向进行，该思想在A.R.Oganov 等开发的晶体结构预测软件USPEX（Universal Structure Predictor：Evolutionary Xtallography）中实现［64］。USPEX通过指纹函数描述晶体结构：

式中，Z为原子相对质量；Rij为两原子间距离，nm；V为单位晶胞体积，nm3；N 为单位晶胞内所含有的原子数，R是一个可变参数。不同结构之间的相似性可由指纹函数空间定义的距离判断：

USPEX 软件基于第一性原理和进化算法，仅需提供材料的化学组分就可以得到能量、体积、硬度、介电常数、能带宽度和磁矩等性质，已经预测了很多无机晶体结构［65-68］，但用于含能晶体预测的报道较少。王洪波［69］用该方法预测了含氮高能量高密度材料的晶体结构，预言了SiC2N4和Si2CN4的高致密相。

3.2.4 粒子群优化算法

粒子群优化算法是在进化算法的基础上发展起来的一种基于种群搜索策略的全局优化算法，它和模拟退火算法相似，都以随机解为起点，通过迭代手段搜索最优解，在该过程中需要定义适应度来标记和评估尝试解的品质。PSO 不需要进行“交叉”和“变异”操作，而直接通过当前搜索的最值来确定下一步的搜索方向，达到对全局最优解的搜索［70］。该算法中粒子更新的速度和位置由以下公式决定：

在结构搜索领域，一些全局性参数优化方法与考虑对称性的结构产生技术、判断相似结构的几何结构因子的引入、以及提高群体多样性的技术等相结合，有效的克服了前期晶体结构预测中难以保持种群多样性的问题，可以保证对势能面进行有效的探索。后来发展的基于高斯过程回归的机器学习势函数，将群体智能与机器学习方法的结合，能进行大尺寸复杂材料的结构预测［73］。

4 基于数学统计方法预测密度

4.1 定量构效关系模型

通过数学和统计学手段找出含能材料的化学结构与密度等性质之间的量变规律，或得到构效关系的数学方程，可以为密度预测提供理论依据［74-76］，核心问题是采用有效算法建立构效关系模型。本文介绍已经用于含能材料密度预测模型的建立方法：回归分析、人工神经网络和改进的遗传算法。

4.1.1 回归分析

回归分析是处理变量之间相关关系的常用数理统计方法［77］，可通过后者的已知或设定值，去估计或预测前者的均值。在回归分析中，如果有两个或两个以上的自变量，就称为多元回归［78］。来蔚鹏等［79］对30 种芳香系单质炸药的偶极矩、最高占据轨道能量、分子总能量等8 个描述符进行多元线性回归计算，建立QSPR 模型预测密度，训练集和测试集的预测值与实验值的平均误差分别为3.33%和2.94%。

支持向量机（Support Vector Machine，SVM）是根据有限的样本信息，在模型复杂性和学习能力之间寻求最佳折中，以获得良好的推广能力。支持向量回归（Support Vector Regression，SVR）是用于解决回归问题的支持向量机。宗朝霞等［74］基于遗传算法的变量筛选和支持向量机预测呋咱系含能化合物的密度，随机选取85%的待测呋咱系含能化合物作为训练集，其余为测试集，预测结果平均相对误差分别为1.16%和2.12%。S.N.Hwang［80］基于多元线性回归和支持向量机建立QSPR 模型预测了3609 种含能材料的密度，结果与Ammo 的基团贡献法相当。

支持向量机和支持向量回归是目前机器学习领域用得较多的方法，序列最小优化算法（Sequential Minimal Optimization，SMO）的提出有效的解决了支持向量机方法实现复杂、效率低等问题，特别在处理大数据量的实际问题时体现出了较好的精度和运算效率，可用于含能材料设计。

4.1.2 人工神经网络

人工神经网络是人工智能研究一种数学模型［81］，在解决复杂环境系统非线性关系模拟等众多问题上展示出强大的功能，越来越多的应用于炸药性能的研究［82-84］。K.Tatyana 等［85］运用人工智能和模式识别方法，通过结构基团分析与识别，从大量已知结构数据库得到规律，对含能材料的性能进行初步预测。人工神经网络可以正确反映配方组成和温度对装药密度的影响规律，对装药密度进行预测［86］。王国栋等［87-88］结合DFT 与人工神经网络，以平均极化率、平均四极矩和偶极矩为描述符，构建炸药密度的预测模型，预测16 种炸药的密度与文献值相对误差为0.47%～6.6%。

神经网络是一个可以自学习并能够“记忆”样本所含信息的理论方法，得到广泛的应用，但它也存在自身的限制与不足，如易形成局部极小值而得不到全局最优解［89］；训练次数多使得学习效率降低，收敛速率慢；隐节点的选取缺乏理论指导等。针对上述问题，国内外已提出不少有效的改进方法，较常用的有增加动量项、自适应调节学习速率、引入陡度因子等。

4.1.3 改进的遗传算法

遗传算法具有自组织、自适应和学习性等特点，搜索过程以目标函数值为评价标准，不受优化函数连续性的约束。针对遗传算法收敛速度慢的缺点，程娥等［90］提出一种伪并行、贪婪的、最优解保存与自适应参数调整相结合的改进遗传算法，对标准化后的399 个炸药分子结构描述符和爆轰数据进行变量的选择并建立定量“分子结构-爆炸性能”关系（QSDR）模型（式（16）），从结构描述符中筛选出13 个与密度相关性最大的参数，呋咱类化合物密度的预测值相对误差在±2%以内。

伪并行是将群体分成N 个子群，各自按照不同的交叉率和变异率进行优化，最后比较各个子群的最优解得到整个群体的最优解；最优解保存策略是每次从所有子群中找出一个最优个体，再将此个体注入每个子群替代各子群中最差个体，防止随机选择时适应值最好个体未被选中的情况；自适应参数调整即交叉率和变异率通过个体本身适应度和种群整体性能的比较确定；贪婪策略通过对父代、子代4 个个体适应度作比较，选择适应度较大的2 个作为子代个体，确保最佳个体被遗传到下一代而不受破坏。改进的遗传算法与传统遗传算法相比，收敛速度快、寻优能力强，误差率小，且性能稳定，有望在含能化合物研发中发挥更大作用。

4.2 经验公式估算法

目前，国际上由经验公式估算含能材料晶体密度的研究中，最有代表性和被广泛引用的报道来自于M.H.Keshavarz。2007 年Keshavarz［91-92］由元素组成预测部分含能化合物的晶体密度，所用公式如下：

MW为分子摩尔质量，g·mol-1；a、b、c、d、ni是C、H、N、O 原子及官能团的数量；经过拟合得到：

对不同的分子结构、不同的官能团及连接方式，还需要分别对上式进行复杂的修正才能用于密度预测，例如对于分子结构中一个C 原子连接三个硝基的情况，密度计算公式如下：

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

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

用多元线性回归方法从实验数据中寻找可调参数z1～z6，EI和ED是与特殊结构相对应的参数，可增加或降低晶体密度，得到密度计算公式：

将72 种化合物的密度预测结果与前期的经验公式法和基团加和法进行比较，准确度有较大提升。

2013 年，Keshavarz［94-97］在对元素组成与特定分子进行修正的同时，综合了定量构效关系方法，考虑分子间的相互作用，应用多元线性回归模型，对含能化合物的密度进行预测，方法的适用性和准确度有较大提高。D.Frem［98］用式（23）的经验关系方法预测的晶体密度值偏高。

CPG和CNG分别是特殊官能团结构对密度的正、负贡献。

2015 年，Keshavarz［99］改进了上述经验关系法，将密度计算公式修正为：

式中，IMP 和DMP 分别是分子间相互作用对密度的修正，体现了特定分子中H 原子数nH和N 原子数nN的偏差，对不同的结构取不同参数，参数的确定来源于实验数据。对172 种不同种类的含能化合物计算得到的均方根误差为2.16%～11.0%。

2017 年，Keshavarz 发布了EMDB（Energetic Materials Designing Bench）软件［100］，可以计算含能材料的30 多种物理化学性质及爆轰性能。

2018 年，V.D.Ghule［101］课题组考虑了含能分子中含有两个及以上NH2或NH3+基团对密度的贡献，将密度计算公式修正为如下形式：

式中，VC、VH、VN、VO是C、H、N、O 的原子体积，nm3。对119 个含苯环的含能材料密度进行计算并与实验数据对比，其中27 个绝对误差小于0.05 g·cm-3，硝基苯及其衍生物中54 个误差范围在3%以内，12 个误差范围是3%～4%，8 个误差范围是5%～8%，离子盐和共晶中25 个误差在3%以内，13 个误差范围是3%～4%，7个误差在5%～7.4%。

QSPR 预测含能晶体密度时，训练样本越多预测效果越好，而实际情况是很多炸药的实验数据缺乏，影响了模型的精度［102-103］，很难明确给出方程的物理意义和分子间相互作用的模式及化学内涵；经验公式方法虽有直接、快速的特点，但普适性差，且忽略了分子间相互作用对晶体堆积结构的影响，在原理上不够准确。

5 展望

综上所述，每种密度计算方法在拥有诸多优势的同时都存在一定的局限性，如基于分子体积的方法在CHNO 类含能材料中应用比较广泛，密度计算的准确度取决于能否正确的描述分子间和分子内相互作用、空间位阻效应、复杂的环结构以及致爆基，如NO2、NH2、OH、N3、ONO2、NHNO2等对密度的影响。MSEP 优于IED 的计算结果，一般会略微低估密度，基组选择对密度的计算结果影响较小。基于力场的三维晶体堆积方法预测的密度一般会比实验值略高，晶体结构搜索要面对巨大的状态空间和高度复杂的能量曲面，发展高效可靠的晶体结构预测方法一直是亟待解决的问题。经验公式法虽然直接但对新发现的分子或某些特殊结构的适用性需验证，QSPR 预测含能晶体密度的模型精度需要提高。以下几种方案结合量子物理、计算化学、人工智能三者的优势，可显著提高含能材料密度计算的效率和可靠性，实现含能材料的高通量设计：

（1）将群体智能搜索与密度泛函理论的色散修正和机器学习势函数方法结合，进行含能晶体的结构和密度预测：先通过高斯过程回归得到势函数，用于势能面快速搜索，提高效率；采用全空间密网格初筛保证初筛样本的完备性，采用刚性分子堆积技术或晶体配位搜索技术减小自由度，并适用于多种成键类型，提高成功率；通过自由能表面全局最小化搜索得到备选晶体结构后，再用DFT-D 方法对晶体结构和密度精细计算，提高置信度。

（2）定量构效关系与人工智能方法结合，可以方便准确的预测含能晶体密度：建立神经网络模型描述复杂的相关关系，结合数据挖掘，可以进一步提高模型的代表性和预测精度；该方案需要高质量的含能材料结构与性能的实验数据做支持，为分子设计和结构-性能关系的预测提供训练依据。

（3）第一原理和机器学习结合（QMML 方法）：通过高通量集成计算产生大批量的含能分子及复合体系，将相关结构和性质的计算结果形成数据库，用于分子表面静电势修正，基于机器学习动态更新样本结构，一方面通过开发专用基函数和高精度赝势，提高第一原理大规模并行计算能力，进行结构构造和性能预估；另一方面，高通量计算与数据挖掘融合，可能发现隐藏的材料“基因”-从原子排列到相结构再到显微组织形成，最后到材料宏观性能与使用寿命之间的关系。