《含能材料》：含能晶体密度预测的研究进展（中国工程物理研究院）

含能晶体密度预测的研究进展 文章编号:1006-9941(2020)01-0001-12 含能晶体密度预测的研究进展 王丽莉12,熊鹰2,谢炜宇2,牛亮亮2,张朝阳 (1.中国工程物理研究院计算机应用研究所,四川绵阳621999;2.中国工程物理研究院化工材料研究所,四川绵阳621999) 摘要:密度是决定含能材料爆轰性能的重要参数。为评估现有CHON类含能材料密度的计算方法,对等电子密度面法、分子表 面静电势法、基团加和法、晶体堆积法、定量构效关系法、经验公式法等进行分析和归类。结果表明,基于分子体积预测方法的精度 取决于分子间和分子内相互作用对密度影响描述的准确度。其中,准确描述氢键和 van der waals作用充满了挑战性。基于晶体体 积计算密度的核心在于晶体结构的准确预测,结构搜索要面对巨大的状态空间和高度复杂的能量曲面的困难,预测效率是亟待解决 的问题。体积加和法和经验公式法存在无法区分同分异构体和晶型的缺点,且对新发现的具有特殊结构的分子由于缺乏实验数据 难以获得准确的经验参数,计算结果偏差较大。引入人工神经网络、遗传算法以及支持向量机等机器学习算法后,定量构效关系法 在含能化合物性能与结构关系研究中取得很大成就,模型精度进一步提高将为基于材料基因组模式的含能材料设计研发奠定基础, 这也是今后密度预测方法发展的主要方向。 关键词:含能材料;密度预测方法;结构预测 中图分类号:T55;064 文献标志码:A DOl:10.11943/CEM2018269 Moller- Plesset perturbation,MP2)理论、耦合簇(Cou 1引言 pled Cluster with Singles and Doubles, CCSD)FI ie 含能材料密度的理论预测方法可分成以下三大等。二是基于晶体体积计算密度,晶体结构预测 类:一是基于分子体积计算密度,由于分子体积没有明 ( Crystal Structure Prediction,CSP)结合了能量和结构 的精确计算方法与全局结构搜索方法,其中结构与能 确的物理定义,因此存在多种计算方法,如等电子密度 量的计算方法包括密度泛函理论的色散校正( Disper 面( Isosurface of Electron Density,ED)方法通常将分 sion Correction Density Functional Theory, DFT-D) 子摩尔体积定义为0.00|a.u.电子密度等值面轮廓范 分子力场与分子动力学方法( Molecular Dynamics, 围内的体积,分子表面静电势( Molecular Surface MD),半经验分子轨道方法( emi-empirical Molecu E| ectrostatic Potentials,MESP)校正法是在密度求解 ar Orbita,MO)等;结构搜索的代表性算法有模拟退 公式中加入通过实验数据拟合的系数表征分子间相互火法( Simulated Annealin,sA)、随机结构搜索法 作用,基团加和( Group Additivity Methods, GAM)法( Ab Initio Random Structure Searching,A|RSs)、进化 是由原子和官能团体积的贡献估算分子体积等,其中算法( Evolutionary Approach,EA)和粒子群优化算法 计算电子密度的方法多数采用密度泛函理论( Density( Particle Swarm Optimization,PsO)等。三是基于数 Functional Theory,DFT)、二阶微扰( Second- order学和统计方法估算含能材料密度,如通过多元线性回 归( Multiple Linear Regressin,MLR)、人工神经网络 稿日期:2018-09-19;修回日期:2019-05-10 网络出版日期:2019-05-24 ( Artificia| Neural networks,ANN)和改进的遗传算法 基金项目:挑战计划(TZ2018004) ( Genetic Algorithm,GA)等建立快速稳定的定量构效 作者简介:王丽莉(1978-),女,副研究员,主要从事计算材料与材 X a( Quantitative Structure Property Relationship 料信息学研究。e-mail;tjhsunwind@osina.com 通信联系人:张朝阳(1971-),男,研究员,主要从事计算含能材料QSPR)模型,或根据经验公式由元素组成和官能团数 学研究。e-mail:chaoyanggang@caep.Cn 目估算晶体密度等。随着云计算的普及和大数据时代 引用本文:王丽莉,熊鷹,谢炜宇,等,含能晶体密度预测的研究进展田.含能材料,2020,28(1)1-12 WANG Li-li. XIONG Ying, XIE Wei-yu et al. Review of Crystal Density Prediction Methods for Energetic Materials[ll. Chinese Joumal of Energetic Materials CHINESE JOURNAL OF ENERGETIC MATERIALS 含能材料 2020年第28卷第1期(1-12

含能材料 2020 年 第 28 卷 第 1 期 （1-12） 含 能 晶 体 密 度 预 测 的 研 究 进 展 CHINESE JOURNAL OF ENERGETIC MATERIALS 含能晶体密度预测的研究进展 王丽莉 1，2 ，熊 鹰 2 ，谢炜宇 2 ，牛亮亮 2 ，张朝阳 2 （1. 中国工程物理研究院计算机应用研究所，四川 绵阳 621999；2. 中国工程物理研究院化工材料研究所，四川 绵阳 621999） 摘 要： 密度是决定含能材料爆轰性能的重要参数。为评估现有 CHON 类含能材料密度的计算方法，对等电子密度面法、分子表 面静电势法、基团加和法、晶体堆积法、定量构效关系法、经验公式法等进行分析和归类。结果表明，基于分子体积预测方法的精度 取决于分子间和分子内相互作用对密度影响描述的准确度。其中，准确描述氢键和 van der Waals 作用充满了挑战性。基于晶体体 积计算密度的核心在于晶体结构的准确预测，结构搜索要面对巨大的状态空间和高度复杂的能量曲面的困难，预测效率是亟待解决 的问题。体积加和法和经验公式法存在无法区分同分异构体和晶型的缺点，且对新发现的具有特殊结构的分子由于缺乏实验数据 难以获得准确的经验参数，计算结果偏差较大。引入人工神经网络、遗传算法以及支持向量机等机器学习算法后，定量构效关系法 在含能化合物性能与结构关系研究中取得很大成就，模型精度进一步提高将为基于材料基因组模式的含能材料设计研发奠定基础， 这也是今后密度预测方法发展的主要方向。 关键词：含能材料；密度预测方法；结构预测 中图分类号：TJ55；O64 文献标志码：A DOI：10.11943/CJEM2018269 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）理论、耦合簇（Cou⁃ pled Cluster with Singles and Doubles，CCSD）理 论 等 。 二 是 基 于 晶 体 体 积 计 算 密 度 ，晶 体 结 构 预 测 （Crystal Structure Prediction，CSP）结合了能量和结构 的精确计算方法与全局结构搜索方法，其中结构与能 量的计算方法包括密度泛函理论的色散校正（Disper⁃ sion Correction Density Functional Theory，DFT⁃D）， 分子力场与分子动力学方法（Molecular Dynamics， MD），半经验分子轨道方法（Semi⁃empirical Molecu⁃ lar 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）模型，或根据经验公式由元素组成和官能团数 目估算晶体密度等。随着云计算的普及和大数据时代 文章编号：1006⁃9941（2020）01⁃0001⁃12 引用本文：王丽莉,熊鹰,谢炜宇,等 . 含能晶体密度预测的研究进展[J]. 含能材料,2020,28(1):1-12. WANG Li‑li, XIONG Ying, XIE Wei‑yu,et al. Review of Crystal Density Prediction Methods for Energetic Materials[J]. Chinese Journal of Energetic Materials （Hanneng Cailiao）,2020,28(1):1-12. 收稿日期：2018⁃09⁃19；修回日期：2019⁃05⁃10 网络出版日期：2019⁃05⁃24 基金项目：挑战计划（TZ2018004） 作者简介：王丽莉（1978-），女，副研究员，主要从事计算材料与材 料信息学研究。e⁃mail：fjhsunwind@sina.com 通信联系人：张朝阳（1971-），男，研究员，主要从事计算含能材料 学研究。e⁃mail：chaoyangzhang@caep.cn 1

王丽莉,熊鹰,谢炜宇,牛亮亮,张朝阳 的到来,机器学习( Machine Learning,ML)及人工智等用 Politzer方法计算了2,4,6-三硝基甲苯 能( Artificial Intelligence,A)技术将推动含能晶体结(TNT)、环三亚甲基三硝胺(RDX)、环四亚甲基四硝胺 构、性质以及感度之间关系的研究,并为开发新型钝感(HMX)等多种含能分子的结构和密度,都与实验数据 高能炸药提供依据。 吻合较好。2013年,B.M.Rice用 Politzer方法对前 期工作进行改进0,并将计算结果分组与实验数据比 2基于分子体积计算密度 较,其中38个中性分子晶体理论密度的平均绝对误差 2.1等电子密度面法 从0.05g·cm减小至0.035g·cm3,48个离子化合物 等电子密度面法计算密度的主要步骤如下:优化 理论密度的平均绝对误差从0.088g·cm3减小至 含能分子确保其处于能量最低点,然后计算振动光谱 0.045g 张君君等0用 Politzer方法计算 和相应热力学信息评估结构稳定性,用 Monte-Carlo N-1,4,6三硝基六氢咪唑(TNNA)和[4,5-d]咪 方法计算每个分子的平均摩尔体积Vam,理论密度唑2(1H)亚硝胺(DNNA晶体密度分别为1.9gcm3 用分子摩尔质量M与vam之比求得口。该方法对分和1.838cm,与实验密度189gcm2和182gcm2 子形状和体积的描述依赖电子密度值的选择,对于处吻合较好。C.K.Kim1将密度计算公式修正为如下 于凝聚态的分子,考虑到分子间压积会使范德华表面形式 有一定穿透,也有用0.002au.或0.0015a.u.的电子 Crystal=1. +0.0189(m)+0.145 AREA 密度等值面作为范德华表面4。 密度泛函理论的B3YP方法和631G“基组是等Ⅱ=∑v(小)- 电子密度面法计算分子晶体密度的主要方法,准确 性优于半经验MO方法間,经济性优于MP2和 式中,AREA是范德华分子表面积,nm2;∏是平均偏差 CCSD(T)等高水平电子密度方法。对CHNO型分V()是静电势,v为分子表面上平均电势,kmor 子晶体,该方法计算中性分子的理论密度与实验值的该方法计算密度的平均绝对误差是0.039g·cm。 平均误差和均方根误差分别为1%和3.7%,分子中含A. Nirwan等将221种炸药分成7类,对比了Lee 氮量超过50%或含氮环的误差为-1.8%和3.4%,离Rice、Kim、 Politzer等MESP方法计算密度的有效性 子晶体的误差为-5.1%和66%,通过对体积公式进行其中Rice和 Politzer修正更接近实验值,多数误差小 修正,离子晶体平均误差和均方根误差减小至-1.1%于3%,硝酸脂和包含笼状或刚性环结构的误差较大 和4.8%10 分子表面静电势校正方法是日前计算含能材料密 等电子密度面法没有考虑分子间空隙及各种复杂度的主要方法,计算精度优于等电子密度面法和基团 的相互作用,用于密度这种宏观性质计算,原理受质加和法,其可靠性己经被证明。 疑。由于该方法对极性分子的计算误差偏大,后期根23基团加和法 据经验进行参数的计算拟合,对不同的分子结构进行 基团加和法用组成分子的原子和官能团的标准体 体积修正,精确度有所提高。 积之和估算分子体积,再用分子摩尔质量求密度。 2.2分子表面静电势校正法 H.L. Ammon等{21从1155个晶体结构数据中提取 P. Politzer等{-1提出修正分子表面静电势参数了78种原子和官能团的体积,建立标准体积数据库 法,在密度求解公式中分别考虑了与分子表面正、负电用于计算 CHNOF类化合物的晶体密度,后来又扩充 荷分布情况相关的分子表面静电势平衡系数v和总方和升级该数据库223,用于更多有机物的密度预测。 差a2,计算公式见式(1) D. Mathieu2以 Ammon数据库为基础,用基团加和法 (1)预测密度与实验数据的平均误差接近2%,只有部分 密堆积结构的估算密度误差达到6%~8%。S.Beau 式中,M为分子摩尔质量,gmo;a、B、y是衰减系camp252用带电基团体积加和法估算离子盐和水合 数,由计算拟合得到。该方法预测密度与实验值误差物的晶体密度,平均误差小于2.5%,通过修正氢键和 般小于0.05gcm3,被广泛应用于CHNO型化合环状结构等对分子体积的贡献,42880个晶体估算 物密度的预测。如H. Singh4、P.Ma1和Y.J.Shu密度的平均误差为2%。 Chinese Joumal of Energetic Materials. Vol28, No 1, 2020(1-12) 含能材料 www.energetic-materials.orgcn

Chinese Journal of Energetic Materials，Vol.28, No.1, 2020（1-12） 含能材料 www.energetic-materials.org.cn 王丽莉，熊鹰，谢炜宇，牛亮亮，张朝阳 的到来，机器学习（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 和总方 差 σ2 total ，计算公式见式（1）： ρ = α ( M V m(0.001) ) + β (vσ2 total) + γ （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］将密度计算公式修正为如下 形式： ρ crystal = 1.06( M AREA ) + 0.0189(Π) + 0.145 （2） Π = 1 n∑i = 1 n | V (ri) - V | ˉ S （3） 式中，AREA 是范德华分子表面积，nm2 ；Π 是平均偏差， V (ri) 是静电势，Vˉ S 为分子表面上平均电势，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. Beau⁃ camp［25-26］ 用带电基团体积加和法估算离子盐和水合 物的晶体密度，平均误差小于 2.5%，通过修正氢键和 环状结构等对分子体积的贡献［27］ ，42880 个晶体估算 密度的平均误差为 2%。 2

含能晶体密度预测的研究进展 基团加和法中体积的贡献值是从大量实验数据中 DFT-D是计算含能材料晶体密度最适合的DFT 总结得出,缺乏实验数据难以获得准确的经验参数。方法,计算精度由泛函和色散校正模型的优良性、两者 不仅如此,该方法不考虑分子构型和分子间相互作用,对色散作用描述的互补性两个因素决定,误差明显小 无法区分同分异构体的密度差异,忽略晶型及温度对于传统的DFT方法。DFTD的可行性与精度经过了 密度的影响,对非电中性含能化合物,含硅、硼等元素系统和全面的验证,但计算量随体系的增大呈三次方 的化合物,或部分密堆积结构含能晶体的计算偏差较或二次方增长,比较耗时。采用对称适应微扰理论 大,估算范围有一定限制 Symmetry Adapted Perturbation Theory, SAPT)JE #H 互作用能直接分解为静电相互作用、诱导能和色散能 3基于晶体结构计算密度 对称性微扰理论定义的交换能,可以计算更大规模的 原子模型。R. Podeswa.用基于DFT单体描述的 3.1量化方法优化晶体结构 从理论上预测含能晶体的结构,可以得到密度。完整势能面 SAPT准确预测了RDX的晶体结构,计算了FOX-7的 晶体结构的求解通常需要高精度量化方法进行结构优311分子动力学结合第一原理 化来实现局域最小值精修和能量稳定性计算。 提高DFT方法计算规模和效率的方法之一是与 3.1.1密度泛函理论的色散校正 分子力场和分子动力学方法结合。MD方法可研究单 密度泛函理论对分子间弱相互作用的描述不完全质炸药及其复合体系的结构、晶体生长机制叫、能量、 准确,为了能够精确的计算含能分子晶体的物性,需要热稳定性、堆积密度、爆轰性能、溶剂对晶体的影响、力 加入分子间的色散校正。非局域性色散校正、Lon 学性能及其温度效应,并与感度关联研究其安全性能, don色散校正、球原子模型等是对DFT的有效修正。计算结果与实验数据越来越接近m 其中DFTD方法在量化程序中得到广泛支持,通常用 分子动力学结合第一原理计算在晶体密度求解方 零阻尼函数(式(5)来调节色散校正在近程、中程距面取得了极大的进展。P. Srinivasan1用B3LYP方法 离时的行为,避免近距离时校正能(式(4))较大导致和6-311C基组优化2,4-二硝基苯甲酸(DNBA)的分 的重复问题2。 子结构,再用 MOLPAK( MOLecular PAcKing)建立空 226 25nRaamp (Rau) (4)间堆积模型,进行UMD力场晶格能最小化“,晶体结 构预测结果与X-ray衍射实验结果的偏差很小,密度的 1+e (5)计算结果与实验值误差为298%。D. N. Kaushik4用 第一原理处理有机晶体中的单分子和短程二体相互作 其中R代表AB原子间的距离,sn是截断半径,m;C用,用极化分子力场处理多体和长程相互作用,建立了 是原子间色散校正系数;sn和s是刻度因子;γ是预设基于片段的杂化多体相互作用模型,晶胞参数计算结 常数。随原子间距离减小,f逐渐为0,故称零阻尼。 果优于色散修正B3LYP-D*和 AMOEBA极化力场方 经过DFTD校正后,静电主导的弱相互作用,如法。.5.L1在B3YP方法和6-31++G*基组计算分 氢键、卤键等问题,计算精度有较大改进,已用于计算子结构的基础上,用 Dreiding力场和模拟退火法预测 5,5′-联四唑-1,1′-二氧二羟铵盐(TKX-50)、1,3,5 了2,4,6,8-四硝基-1,3,5,7-四氮杂环丁烷 氨基2,4,6-三硝基苯(TATB)、1,1-二氨基-2,2-二硝( TNTAZC)和2,3,5,6,7,8-六硝基-1,4-二氮杂环丁 基乙烯(FOX-7)、2,6-二氨基-3,5-二硝基吡嗪-1-氧烷( HNDAZC)的晶体密度分别为2.156g·cm3和 (LLM-105)、BHMX和六硝基六氮杂异戊兹烷1.975gcm-3。力场选择对晶体堆积和密度计算非常 (CL-20)等晶体密度,晶胞参数计算结果与实验的误关键, Dreiding力场对84%的不含芳香环结构的密度 差通常在2%以内。L. Zhang等采用PBE交换关联预测值与实验值的偏差在5%以内,PCF力场对83% 势计算HMX的a、B、δ和y相的晶体结构,体积误差分别的含芳香环结构的密度预测值与实验值的偏差在4% 为-2.36%,-1.00%,+0.12%和-0.3%,采用vdWD3以内,广义 Amber力场计算FOX-7、RDX、CL-20 方法计算了53种炸药晶体结构,体积误差<±4%。HMX等的密度比实验值低5%~8%41,CVF力场和 H.Lin3、H.C.S.Chan13、P. Panin1也采用 DFT-D DMACRYS力场也用于含芳香环结构的含能晶体预 对含能晶体或共晶进行结构和密度预测 测,密度误差小于7%4-47 CHINESE JOURNAL OF ENERGETIC MATERIALS 含能材料 2020年第28卷第1期(1-12

CHINESE JOURNAL OF ENERGETIC MATERIALS 含能材料 2020 年 第 28 卷 第 1 期 （1-12） 含 能 晶 体 密 度 预 测 的 研 究 进 展 基团加和法中体积的贡献值是从大量实验数据中 总结得出，缺乏实验数据难以获得准确的经验参数。 不仅如此，该方法不考虑分子构型和分子间相互作用， 无法区分同分异构体的密度差异，忽略晶型及温度对 密度的影响，对非电中性含能化合物，含硅、硼等元素 的化合物，或部分密堆积结构含能晶体的计算偏差较 大，估算范围有一定限制。 3 基于晶体结构计算密度 3.1 量化方法优化晶体结构 从理论上预测含能晶体的结构，可以得到密度。 晶体结构的求解通常需要高精度量化方法进行结构优 化来实现局域最小值精修和能量稳定性计算。 3.1.1 密度泛函理论的色散校正 密度泛函理论对分子间弱相互作用的描述不完全 准确，为了能够精确的计算含能分子晶体的物性，需要 加入分子间的色散校正。非局域性色散校正、Lon⁃ don 色散校正、球原子模型等是对 DFT 的有效修正。 其中 DFT⁃D 方法在量化程序中得到广泛支持，通常用 零阻尼函数（式（5））来调节色散校正在近程、中程距 离时的行为，避免近距离时校正能（式（4））较大导致 的重复问题［28］ 。 E DFT - D3 disp = - 1 2 ∑A ≠ B n ∑= 6.8 sn C AB n R n AB fdamp，n (RAB ) （4） fdamp，n (RAB ) = é ë ê1 + e -γ (RAB /sr，nR ) AB 0 - 1 ù û ú -1 （5） 其中 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］ 。 3

王丽莉,熊鹰,谢炜宇,牛亮亮,张朝阳 分子动力学的核心在于求解正则运动方程,分子发展了许多晶体结构预测的方法,为研究含能晶体的 内和分子间的作用力都包含在分子力场中,因此计算结构及其性质注入活力。2016年,多个单位参加了对 的准确性依赖所使用的力场4-。如果缺乏完全适有机晶体结构的第六次盲测2,分子结构均来自第 用的力场参数,无法准确描述分子间的相互作用,致使原理计算和剑桥结构数据库,晶体结构通过随机搜索 晶体结构预测存在较大偏差,要提高晶体密度计算的法、遗传算法或模拟退火法产生。 准确性需自研发经验势。 3.2.1模拟退火法 313半经验方法 模拟退火法是S. Kirkpatrick等于1983年基于 半经验方法通过忽略分子积分计算或引进经验参金属退火机理而建立的全局最优化方法,通过对给定 量对 Hartree-Fock方程的近似求解,如电荷自恰密度初始状态进行模拟加温,生成一系列关联的新解并评 泛函紧束缚( Self-consistent Charge Density Function-估与初始解的能量差,再根据 Metropolis接受准则进 al Tight-binding, SCC-DFTB)方法是在DFT基础上采行结构筛选,得到一系列给定温度下允许存在新 用非正交基的非正交化方式来对紧束缚( Tight Bind-解。例如在温度T时,由当前状态i产生新状态能量 nTB)理论做修正,计算速度比传统DFT快两三个分别为E和E,若E>E则接受新状态为当前状态,否 数量级,核心是对哈密顿矩阵H的求解和对电子波函则,以概率P来接受状态,其中k为 Boltzmann常数。 数的计算,能量表达式如下所示 p=exp[ -(E-E)/k] Em=∑{的HP)+∑E (6) 在模拟退火的开始阶段通常以升高温度来实现初 式中,Oc表示电子占据数,i表示具体的价电子,EA始解在能量曲面上的上坡游走,克服势垒走向周围的 表示原子核A和B之间的排斥能,k。电子波函数主 局域最小值;这些走离初始解位置的新解通常都处于 要通过基函数的线性组合来求解: 能量面中的坡上,需要进一步通过降温处理让它们再 v(x,t)=∑va,t)(x) (7)次在能量曲面上进行下陂游走,走向周围的局域态 升温和降温时,温度变化遵循如下关系式 式中,l是电子的空间坐标,a=x,y,z表示空间维度 Tnew= Told x(1.0+T) 大大降低了哈密顿矩阵的求解难度,可用于计算晶体m=mx(0-7) (10) 结构和密度0,但该方法没有针对硼元素的经验参式中,Tm和Ta分别是变温前后的温度,K;T和T分 数,无法模拟极端高温高压条件下的物理化学变化。别是自定义的加热因子和冷却因子。 Holden用半经验AM1( Austin Model number1) Yin153和G.Z.Zhao5分别用模拟退火技术预 方法对堆积结构进行优化,预测晶胞参数。 测了1′-氨基-1H-1,5′双四唑(ATT)及其溴氰衍生物 半经验方法虽然可将计算规模增大,但还不能有和2,4,6三硝基-1,3,5-三氮1,3,5-三氧化六元氮杂 效以及长时间模拟原子数量非常庞大的系统或者结构环( TNTATO)的晶体结构。刘英哲等”基于B3LY 更加灵活的体系,比如复杂分子表面势能以及弱键能方法和6-31G(d)基组建立了适用于全氮结构的力场, 等;计算的可靠性与参数获得的途径密切相关,参数拟用 Monte-Carlo方法进行空间堆积,预测晶体结构并计 合大多依赖于DFT计算;由于采用了一系列近似,应算N4、N4、Na、N0、N12五种笼型全氮分子的晶体密度分 用于含有强吸电子基团(如硝基)的体系时,由于电荷别是1.81,2.08,2.47,2.46,2.57gcm。.F 密度不对称,计算误差会大大增加 Moxnes用 Polymorph Predictor方法和 COMPASS 3.2势能面搜索预测晶体结构 力场计算晶体密度并与实验数据比较(见表1),对大 晶体结构预测是在固定化学组分情况下确定最有多数空间群,模拟退火法高估了密度,平均绝对误差 可能的原子排列方式和晶胞参数,这一过程中用到自0.07~0.08gcm3。 然界中占据统治地位的能量最低原理,即在固定组分 模拟退火法是早期含能晶体结构预测常用方法, 的状态子空间中搜索处于能量曲面极小值的排列方这种基于初始结构演化的跃迁势垒方法对人为选择的 式。从理论上预测晶体结构一直是化学、物理和材料初始结构具有依赖性,考虑空间群不够宽泛,难以克服 科学中的一大难题,随着理论方法和计算技术的完善,退火过程的高能量势垒,需要进行多次迭代直到所得 Chinese Joumal of Energetic Materials. Vol28, No 1, 2020(1-12) 含能材料 www.energetic-materials.orgcn

Chinese Journal of Energetic Materials，Vol.28, No.1, 2020（1-12） 含能材料 www.energetic-materials.org.cn 王丽莉，熊鹰，谢炜宇，牛亮亮，张朝阳 分子动力学的核心在于求解正则运动方程，分子 内和分子间的作用力都包含在分子力场中，因此计算 的准确性依赖所使用的力场［48-49］。如果缺乏完全适 用的力场参数，无法准确描述分子间的相互作用，致使 晶体结构预测存在较大偏差，要提高晶体密度计算的 准确性需自研发经验势。 3.1.3 半经验方法 半经验方法通过忽略分子积分计算或引进经验参 量对 Hartree⁃Fock 方程的近似求解，如电荷自恰密度 泛函紧束缚（Self⁃consistent Charge Density Function⁃ al Tight⁃binding，SCC⁃DFTB）方法是在 DFT 基础上采 用非正交基的非正交化方式来对紧束缚（Tight Bind⁃ ing，TB）理论做修正，计算速度比传统 DFT 快两三个 数量级，核心是对哈密顿矩阵 H 的求解和对电子波函 数 λ的计算，能量表达式如下所示： EDFTB =∑i occ ψi | H | 0 ψi + ∑A ≠ B atoms E AB rep （6） 式中，occ 表示电子占据数，i 表示具体的价电子，E AB rep 表示原子核 A 和 B 之间的排斥能，kJ。电子波函数主 要通过基函数的线性组合来求解： Ψi (x，t ) =∑la Ψi (la，t )Φla (x ) （7） 式中，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 常数。 p = exp[ - (Ej - Ei) /kt] （8） 在模拟退火的开始阶段通常以升高温度来实现初 始解在能量曲面上的上坡游走，克服势垒走向周围的 局域最小值；这些走离初始解位置的新解通常都处于 能量面中的坡上，需要进一步通过降温处理让它们再 次在能量曲面上进行下陂游走，走向周围的局域态。 升温和降温时，温度变化遵循如下关系式： Tnew = Tol d × (1.0 + Th ) （9） Tnew = Tol d × (1.0 - Tc ) （10） 式中，Tnew 和 Tol d 分别是变温前后的温度，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 。 模拟退火法是早期含能晶体结构预测常用方法， 这种基于初始结构演化的跃迁势垒方法对人为选择的 初始结构具有依赖性，考虑空间群不够宽泛，难以克服 退火过程的高能量势垒，需要进行多次迭代直到所得 4

含能晶体密度预测的研究进展 结构满足所设置的收敛标准,在计算效率上有很大3.2.3进化算法 缺陷 进化算法是通过不断地产生候选晶体结构种群, 并经过生物进化理论中的遗传繁衍、基因变异、基因重 表1 Polymorph Predictor方法计算的密度与实验数据对比组和自然选择等基本法则不断更新和优化这一种群 Table 1 Comparison of the calculated densities by poly. 直至搜索到满足要求的晶体结构。在这一过程中,需 morph Predictor methods and experimental data 要设定一个更新法则来决定如何进行繁衍、变异和重 MRE /g·cm 组,同时还需要一个适合的选择法则来决定个体的存 RDX P2,c 活几率、变异几率和杂交几率,以保证个体进化能够朝 P21c10 着优化的方向进行,该思想在A.R. Oganov等开发的 NTO P21/c 晶体结构预测软件 USPEX( Universal Structure Pre DNAM dictor: Evolutionary Xtallography)中实现。 USPEX E-CL-202.01[5】P2 通过指纹函数描述晶体结构 .8951P21212112 X1.8759P2 P21/c )=∑∑22a(-R) (12) a-HMX1.96151C2/c Note:1) Pesn is experimental densities). is calculated space groups of式中,z为原子相对质量;Rn为两原子间距离,nm;V为 minimum in total energy, and G is calculated space groups of mini- 单位晶胞体积,nm3;N为单位晶胞内所含有的原子 when space group is G, and MRE, is deviation of calculated den.数,R是一个可变参数。不同结构之间的相似性可由 sities when space group is Ga respectively 指纹函数空间定义的距离判断 3.2.2随机搜索法 FPFP da=0.51 (13) 随机搜索方法通过在状态空间尽可能多地进行无 FPⅢFP 区别随机采样,产生大量初始结构X,参照(11)式在 其周围以一定的步长随机搜索可行解 (R) f(R)dR (14) X,=X,+(b-a).rand( (11) USPEX软件基于第一性原理和进化算法,仅需提 式中,t为当前迭代次数,rand()为随机数函数,a和b供材料的化学组分就可以得到能量、体积、硬度、介电 为随机数的取值范围;然后进行精细结构优化使它们常数、能带宽度和磁矩等性质,已经预测了很多无机晶 走向邻近的能量极小值,经过不断迭代从而完成晶体体结构65,但用于含能晶体预测的报道较少。王洪 结构的预测61。 波用该方法预测了含氮高能量高密度材料的晶体 第一原理随机搜索方法将DFT-D和ARSS结合,在结构,预言了SiC2N4和SiCN4的高致密相。 晶体结构预测中显现出了强大的能力。如M.zika23.2.4粒子群优化算法 用该方法成功的预测了有机分子m-氨基苯甲酸 粒子群优化算法是在进化算法的基础上发展起来 (mABA)的晶体结构,与Xray衍射和固态核磁共振的一种基于种群搜索策略的全局优化算法,它和模拟 的低温实验结果吻合,并能区分不同类型氢键的作用。退火算法相似,都以随机解为起点,通过迭代手段搜索 T.J.Linω用ARSs方法发现了甲醇新的δ相候选结最优解,在该过程中需要定义适应度来标记和评估尝 构,通过DFT-D计算确定了相变压力和CH…O氢键试解的品质。PSO不需要进行“交叉”和“变异”操作 在高压下晶体结构的稳定性中起主导作用。 而直接通过当前搜索的最值来确定下一步的搜索方 随机搜索方法的优点是对目标函数无特殊要求,向,达到对全局最优解的搜索0。该算法中粒子更新 通用性强。但计算需要足够大的样本方可获得较精确的速度v和位置x由以下公式决定 的结构。依据品格能的计算得到品体结构之间的能量=wv+(pnm-刘)+c(8m,-对)x”(15) 差很小,随机搜索方法仅按照能量进行筛选很容易漏 =x:+y 掉其他热力学稳定结构。在搜索进程中增加自适应功式中,v表示粒子当前速度向量,ms;x表示粒子当 能,使算法能够根据当前解自动计算下次迭代搜索的前位置向量;w为惯性权重,用于调整历史速度对当前 方向和步长,可有效改进搜索能力和精度。 速度的影响,平衡局部搜索和全局搜索能力;P=m表 CHINESE JOURNAL OF ENERGETIC MATERIALS 含能材料 2020年第28卷第1期(1-12

CHINESE JOURNAL OF ENERGETIC MATERIALS 含能材料 2020 年 第 28 卷 第 1 期 （1-12） 含 能 晶 体 密 度 预 测 的 研 究 进 展 结构满足所设置的收敛标准，在计算效率上有很大 缺陷。 3.2.2 随机搜索法 随机搜索方法通过在状态空间尽可能多地进行无 区别随机采样，产生大量初始结构 Xst ，参照（11）式在 其周围以一定的步长随机搜索可行解： Xt = Xst + (b - a) ⋅ ran d( ) + a （11） 式中，t 为当前迭代次数，ran d( )为随机数函数，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 Pre⁃ dictor：Evolutionary Xtallography）中实现［64］ 。USPEX 通过指纹函数描述晶体结构： f (R) =∑i ∑ j ≠ i ZiZj 4πR 2 ij V N δ (R - Rij) （12） 式中，Z 为原子相对质量；Rij为两原子间距离，nm；V 为 单位晶胞体积，nm3 ；N 为单位晶胞内所含有的原子 数，R 是一个可变参数。不同结构之间的相似性可由 指纹函数空间定义的距离判断： dij = 0.5 (1 - FPi FPi  FPi FPi ) （13） FPi(R) = 1 D ∫iD (i + 1)D f n (R) d R （14） USPEX 软件基于第一性原理和进化算法，仅需提 供材料的化学组分就可以得到能量、体积、硬度、介电 常数、能带宽度和磁矩等性质，已经预测了很多无机晶 体结构［65-68］ ，但用于含能晶体预测的报道较少。王洪 波［69］用该方法预测了含氮高能量高密度材料的晶体 结构，预言了 SiC2N4和 Si2CN4的高致密相。 3.2.4 粒子群优化算法 粒子群优化算法是在进化算法的基础上发展起来 的一种基于种群搜索策略的全局优化算法，它和模拟 退火算法相似，都以随机解为起点，通过迭代手段搜索 最优解，在该过程中需要定义适应度来标记和评估尝 试解的品质。PSO 不需要进行“交叉”和“变异”操作， 而直接通过当前搜索的最值来确定下一步的搜索方 向，达到对全局最优解的搜索［70］ 。该算法中粒子更新 的速度 v t + 1 i，j 和位置 x t + 1 ij 由以下公式决定： v t + 1 i，j = wv t i，j + c1r1 (pbest t i，j - x t ij) + c2r2 (gbest t i，j - x t ij) x t + 1 ij = x t ij + v t + 1 i，j （15） 式中，v t i,j 表示粒子当前速度向量，m·s -1 ；x t ij 表示粒子当 前位置向量；w 为惯性权重，用于调整历史速度对当前 速度的影响，平衡局部搜索和全局搜索能力；pbest t i,j 表 表 1 Polymorph Predictor 方法计算的密度与实验数据对比 Table 1 Comparison of the calculated densities by Poly⁃ morph Predictor methods and experimental data model RDX TNT NTO DNAM ε⁃CL⁃20 DADNE β⁃HMX α⁃HMX ρexp / g·cm-3 1） 1.82［59］ 1.65［59］ 1.91［59］ 2.00［60］ 2.01［59］ 1.89［59］ 1.87［59］ 1.96［59］ Ge 2） P21/c Pbca P21/c Pnma P21 P212121 P21/c C2/c MREe / %3） 3 14 8 4 6 12 0.4 1 Gd 2） P21 P21/c Pbca P212121 P1 Pbca P21/c C2/c MREd / %3） 1 10 5 0.5 2 1 0.4 1 Note： 1）ρexp is experimental densities. 2）Ge is calculated space groups of minimum in total energy，and Gd is calculated space groups of mini⁃ mum deviation in densities. 3）MREe is deviation of calculated densi⁃ ties when space group is Ge，and MREd is deviation of calculated den⁃ sities when space group is Gd respectively. 5

王丽莉,熊鹰,谢炜宇,牛亮亮,张朝阳 示个体最优位置,8=表示全局最优位置;c表示加速向量机建立QSPR模型预测了3609种含能材料的密 步长,c1为控制粒子飞向个体最优位置的步长,c2为控度,结果与Ammo的基团贡献法相当 制粒子飞向全局最优位置的步长,r为独立的随机函 支持向量机和支持向量回归是目前机器学习领域 数,t为代数。这种算法在精度和收敛速度上都有明显用得较多的方法,序列最小优化算法( Sequential Min- 优点,晶体结构搜索软件 CALYPSO( Crystal! tructure imal Optimization,SMO)的提出有效的解决了支持 Analysis Particle Swarm Optimization)实现了该方向量机方法实现复杂效率低等问题,特别在处理大数 法{-21,其高效可靠性经过大量系统的测试。 据量的实际问题时体现出了较好的精度和运算效率 在结构搜索领域,一些全局性参数优化方法与考可用于含能材料设计。 虑对称性的结构产生技术、判断相似结构的几何结构4.1.2人工神经网络 因子的引入、以及提高群体多样性的技术等相结合,有 人工神经网络是人工智能研究一种数学模型。1, 效的克服了前期晶体结构预测中难以保持种群多样性在解决复杂环境系统非线性关系模拟等众多问题上展 的问题,可以保证对势能面进行有效的探索。后来发示出强大的功能,越来越多的应用于炸药性能的研 展的基于高斯过程回归的机器学习势函数,将群体智宽20。K. Tatyana等运用人工智能和模式识别方 能与机器学习方法的结合,能进行大尺寸复杂材料的法,通过结构基团分析与识别,从大量已知结构数据库 结构预测切。 得到规律,对含能材料的性能进行初步预测。人工神 经网络可以正确反映配方组成和温度对装药密度的影 基于数学统计方法预测密度 响规律,对装药密度进行预测。。王国栋等结合 DFT与人工神经网络,以平均极化率、平均四极矩和偶 4.1定量构效关系模型 极矩为描述符,构建炸药密度的预测模型,预测16种 通过数学和统计学手段找出含能材料的化学结构 与密度等性质之间的量变规律,或得到构效关系的数炸药的密度与文献值相对误差为047%-6.6%。 学方程,可以为密度预测提供理论依据4-61,核心问 神经网络是一个可以自学习并能够“记忆”样本所 题是采用有效算法建立构效关系模型。本文介绍已经 含信息的理论方法,得到广泛的应用,但它也存在自身 的限制与不足,如易形成局部极小值而得不到全局最 用于含能材料密度预测模型的建立方法:回归分析、人 优解;训练次数多使得学习效率降低,收敛速率慢 工神经网络和改进的遗传算法 隐节点的选取缺乏理论指导等。针对上述问题,国内 4.1.1回归分析 回归分析是处理变量之间相关关系的常用数理统外已提出不少有效的改进方法,较常用的有增加动量 计方法,可通过后者的已知或设定值,去估计或预 项、自适应调节学习速率、引入陡度因子等 测前者的均值。在回归分析中,如果有两个或两个以 4.1.3改进的遗传算法 上的自变量,就称为多元回归。来蔚鹏等切对 遗传算法具有自组织、自适应和学习性等特点,搜 索过程以目标函数值为评价标准,不受优化函数连续 30种芳香系单质炸药的偶极矩、最高占据轨道能量 分子总能量等8个描述符进行多元线性回归计算,建性的约束。针对遗传算法收敛速度慢的缺点,程娥 立QsPR模型预测密度,训练集和测试集的预测值与 等提出一种伪并行、贪婪的、最优解保存与自适应 参数调整相结合的改进遗传算法,对标准化后的 实验值的平均误差分别为3.33%和294%。 支持向量机( Support Vector Machine,SWM)是399个炸药分子结构描述符和爆轰数据进行变量的选 择并建立定量“分子结构-爆炸性能”关系(QSDR)模 根据有限的样本信息,在模型复杂性和学习能力之间 型(式(16)),从结构描述符中筛选出13个与密度相关 寻求最佳折中,以获得良好的推广能力。支持向量回 归( Support Vector Regression,SVR)是用于解决回归性最大的参数,呋咱类化合物密度的预测值相对误差 问题的支持向量机。宗朝霞等基于遗传算法的变在±2%以内 量筛选和支持向量机预测呋咱系含能化合物的密度,y=0.3873x,-0.972x0+0.4696X0-0.0711X 0.6778x1a2+0.7236X16+0.2086X194+ 随机选取85%的待测呋咱系含能化合物作为训练集 (16) 0.2248x216-0.6644X22+0.1238X239 其余为测试集,预测结果平均相对误差分别为1.16% 0.3443X2x-0.8852X130+0.1952X38 和212%。S.N. Hwang基于多元线性回归和支持 伪并行是将群体分成N个子群,各自按照不同的 Chinese Joumal of Energetic Materials. Vol28, No 1, 2020(1-12) 含能材料 www.energetic-materials.orgcn

Chinese Journal of Energetic Materials，Vol.28, No.1, 2020（1-12） 含能材料 www.energetic-materials.org.cn 王丽莉，熊鹰，谢炜宇，牛亮亮，张朝阳 示个体最优位置，g best t i,j 表示全局最优位置；c 表示加速 步长，c 1 为控制粒子飞向个体最优位置的步长，c 2 为控 制粒子飞向全局最优位置的步长，r 为独立的随机函 数，t为代数。这种算法在精度和收敛速度上都有明显 优点，晶体结构搜索软件 CALYPSO（Crystal Structure Analysis Particle Swarm Optimization）实 现 了 该 方 法［71-72］ ，其高效可靠性经过大量系统的测试。 在结构搜索领域，一些全局性参数优化方法与考 虑对称性的结构产生技术、判断相似结构的几何结构 因子的引入、以及提高群体多样性的技术等相结合，有 效的克服了前期晶体结构预测中难以保持种群多样性 的问题，可以保证对势能面进行有效的探索。后来发 展的基于高斯过程回归的机器学习势函数，将群体智 能与机器学习方法的结合，能进行大尺寸复杂材料的 结构预测［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 Min⁃ imal 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% 以内。 y=0.3873x 9-0.9772x 40+0.4696x 60-0.0711x 107- 0.6778x 182+0.7236x 186+0.2086x 194+ 0.2248x 216-0.6644x 223+0.1238x 239- 0.3443x 270-0.8852x 330+0.1952x 338 （16） 伪并行是将群体分成 N 个子群，各自按照不同的 6

含能晶体密度预测的研究进展 交叉率和变异率进行优化,最后比较各个子群的最优式法和基团加和法进行比较,准确度有较大提升 解得到整个群体的最优解;最优解保存策略是每次从 2013年, Keshavarz9在对元素组成与特定分 所有子群中找出一个最优个体,再将此个体注入每个子进行修正的同时,综合了定量构效关系方法,考虑分 子群替代各子群中最差个体,防止随机选择时适应值子间的相互作用,应用多元线性回归模型,对含能化合 最好个体未被选中的情况;自适应参数调整即交叉率物的密度进行预测,方法的适用性和准确度有较大提 和变异率通过个体本身适应度和种群整体性能的比高。D.Frem9用式(23)的经验关系方法预测的晶体 较确定;贪婪策略通过对父代、子代4个个体适应度密度值偏高。 作比较,选择适应度较大的2个作为子代个体,确保 29.20a+1.515b+53.06c+61.30d =-1.609+ 最佳个体被遗传到下一代而不受破坏。改进的遗传 算法与传统遗传算法相比,收敛速度快、寻优能力强 0.0703Cpc-0.0751CNc (23 误差率小,且性能稳定,有望在含能化合物研发中发 Cc和Cc分别是特殊官能团结构对密度的正、负 挥更大作用。 贡献。 4.2经验公式估算法 2015年, Keshavarz9改进了上述经验关系法,将 目前,国际上由经验公式估算含能材料晶体密度密度计算公式修正为 的研究中,最有代表性和被广泛引用的报道来自于 =1.753-10.238n+9.908nN+0.0992MP M.H. Keshavarz。2007年 Keshavarz.92由元素组成 0.0845DMP (24) 预测部分含能化合物的晶体密度,所用公式如下 式中,MP和DMP分别是分子间相互作用对密度的修 Aa+Bb+Cc+Dd+∑n,(FG) (17)正,体现了特定分子中H原子数n和N原子数n的偏 差,对不同的结构取不同参数,参数的确定来源于实验 Mw为分子摩尔质量,gmor;a、b、c、d、n,是 数据。对172种不同种类的含能化合物计算得到的均 H、N、O原子及官能团的数量;经过拟合得到 方根误差为2.16%~11.0% =(479732a-19.2948b+26.5340C+26.0011d 2017年, Keshavarz发布了EMDB( Energetic Ma 253190n。-0.6358n+11.5413no)/M(18) terials Designing Bench)软件0),可以计算含能材料 对不同的分子结构、不同的官能团及连接方式,还的30多种物理化学性质及爆轰性能。 需要分别对上式进行复杂的修正才能用于密度预测 2018年,V.D.Ghue10课题组考虑了含能分子 例如对于分子结构中一个C原子连接三个硝基的情中含有两个及以上NH2或NH基团对密度的贡献 况,密度计算公式如下 将密度计算公式修正为如下形式 ncH>1.5nNo. Po. corr =3.033-Po (19) +0.00175 (25) av+ bv.+ cv.+ dv ncH.≤0.6no, p=-0.3788+1.2569。(20) 式中,V、V、V、V是C、H、N、O的原子体积,nm3。 其余情况用poc=0.1233+0.83730修正。 对119个含苯环的含能材料密度进行计算并与实验数 2009年, Keshavarz对前期建立的密度计算公据对比,其中27个绝对误差小于0.05g:cm-3,硝基苯 式修正如下 及其衍生物中54个误差范围在3%以内,12个误差范 p=21+2M÷、bx1M=sE,+zED(21)围是3%-4%,8个误差范围是5%~8%,离子盐和共晶 中25个误差在3%以内,13个误差范围是3%~4 用多元线性回归方法从实验数据中寻找可调参数个误差在5%~7.4%。 z1~z,E和E是与特殊结构相对应的参数,可增加或 QSPR预测含能晶体密度时,训练样本越多预测效 降低晶体密度,得到密度计算公式 果越好,而实际情况是很多炸药的实验数据缺乏,影响 p=1.521+6.946 11.53—+20.10 了模型的精度0-103,很难明确给出方程的物理意义和 Mw 分子间相互作用的模式及化学内涵;经验公式方法虽 0.1325E1-0.1559E (22)有直接、快速的特点,但普适性差,且忽略了分子间相 将72种化合物的密度预测结果与前期的经验公互作用对晶体堆积结构的影响,在原理上不够准确。 CHINESE JOURNAL OF ENERGETIC MATERIALS 含能材料 2020年第28卷第1期(1-12

CHINESE JOURNAL OF ENERGETIC MATERIALS 含能材料 2020 年 第 28 卷 第 1 期 （1-12） 含 能 晶 体 密 度 预 测 的 研 究 进 展 交叉率和变异率进行优化，最后比较各个子群的最优 解得到整个群体的最优解；最优解保存策略是每次从 所有子群中找出一个最优个体，再将此个体注入每个 子群替代各子群中最差个体，防止随机选择时适应值 最好个体未被选中的情况；自适应参数调整即交叉率 和变异率通过个体本身适应度和种群整体性能的比 较确定；贪婪策略通过对父代、子代 4 个个体适应度 作比较，选择适应度较大的 2 个作为子代个体，确保 最佳个体被遗传到下一代而不受破坏。改进的遗传 算法与传统遗传算法相比，收敛速度快、寻优能力强， 误差率小，且性能稳定，有望在含能化合物研发中发 挥更大作用。 4.2 经验公式估算法 目前，国际上由经验公式估算含能材料晶体密度 的研究中，最有代表性和被广泛引用的报道来自于 M. H. Keshavarz。2007 年 Keshavarz［91⁃92］ 由元素组成 预测部分含能化合物的晶体密度，所用公式如下： ρ0 = Aa + Bb + Cc + D d +∑ni (FG ) i M W （17） M W 为分子摩尔质量，g·mol-1 ；a、b、c、d、ni 是 C、 H、N、O 原子及官能团的数量；经过拟合得到： ρ0 =(47.9732a - 19.2948b + 26.5340c + 26.0011 d - 25.3190ncoo - 0.6358no + 11.5413n OH )/M W （18） 对不同的分子结构、不同的官能团及连接方式，还 需要分别对上式进行复杂的修正才能用于密度预测， 例如对于分子结构中一个 C 原子连接三个硝基的情 况，密度计算公式如下： nCH2 ≥ 1.5nNO2 ，ρ0，corr = 3.033 - ρ0 （19） nCH2 ≤ 0.6nNO2 ，ρ0，corr = -0.3788 + 1.2569ρ0 （20） 其余情况用 ρ0，corr = 0.1233 + 0.8373ρ0 修正。 2009 年，Keshavarz［93］对前期建立的密度计算公 式修正如下： ρ = z 1 + z 2 a M W + z 3 b M W + z 4 c M W + z 5EI + z 6ED （21） 用多元线性回归方法从实验数据中寻找可调参数 z1~z6，EI和 ED 是与特殊结构相对应的参数，可增加或 降低晶体密度，得到密度计算公式： ρ = 1.521 + 6.946 a M W - 11.53 b M W + 20.10 c M W + 0.1325EI - 0.1559ED （22） 将 72 种化合物的密度预测结果与前期的经验公 式法和基团加和法进行比较，准确度有较大提升。 2013 年，Keshavarz［94-97］在对元素组成与特定分 子进行修正的同时，综合了定量构效关系方法，考虑分 子间的相互作用，应用多元线性回归模型，对含能化合 物的密度进行预测，方法的适用性和准确度有较大提 高。D. Frem［98］ 用式（23）的经验关系方法预测的晶体 密度值偏高。 ρ0 = -1.609 + 29.20a + 1.515b + 53.06c + 61.30 d M W + 0.0703CPG - 0.0751C NG （23） CPG 和 C NG 分别是特殊官能团结构对密度的正、负 贡献。 2015 年，Keshavarz［99］ 改进了上述经验关系法，将 密度计算公式修正为： ρ = 1.753 - 10.238nH + 9.908nN + 0.0992IMP - 0.0845DMP （24） 式中，IMP 和 DMP 分别是分子间相互作用对密度的修 正，体现了特定分子中 H 原子数 nH和 N 原子数 nN的偏 差，对不同的结构取不同参数，参数的确定来源于实验 数据。对 172 种不同种类的含能化合物计算得到的均 方根误差为 2.16%~11.0%。 2017 年，Keshavarz 发布了 EMDB（Energetic Ma⁃ terials Designing Bench）软件［100］ ，可以计算含能材料 的 30 多种物理化学性质及爆轰性能。 2018 年，V. D. Ghule［101］ 课题组考虑了含能分子 中含有两个及以上 NH2 或 NH3 +基团对密度的贡献， 将密度计算公式修正为如下形式： ρ0 = M W aVC + bV H + cV N + d V O + 0.00175 （25） 式中，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］ ，很难明确给出方程的物理意义和 分子间相互作用的模式及化学内涵；经验公式方法虽 有直接、快速的特点，但普适性差，且忽略了分子间相 互作用对晶体堆积结构的影响，在原理上不够准确。 7

王丽莉,熊鹰,谢炜宇,牛亮亮,张朝阳 参考文献 5展望 [1] Bader R F W. Carroll M T, Cheeseman IR, et al. Properties of atoms in molecule: atomic volumes [] Journal of America 综上所述,每种密度计算方法在拥有诸多优势的 Chemical Society,1987,109(26):7968-7979. [2] Rahhi S. Hari J S, Sengupta S K. Computational studies on 同时都存在一定的局限性,如基于分子体积的方法在 2.4-Triazolium-based salts as energetic materials[)].Jouma CHNO类含能材料中应用比较广泛,密度计算的准确 of Chemical Sciences. 2015, 127(6):1099-1107 度取决于能否正确的描述分子间和分子内相互作用、 [3]张驰,陈沫,陈潮,等.稠环类1,2,4,5-四嗪衍生物结构和性能 的理论研究[门]含能材料,2017,25(4):273-281 空间位阻效应、复杂的环结构以及致爆基,如NO2 ZHANG Chi. CHen MO. Chen Chao, et al theoretical NH2、OH、N1、ONO2、NHNO2等对密度的影响。 study on structure and properties of polycyclic derivatives of MSEP优于仼D的计算结果,一般会略微低估密度,基 1.2.4.5-Tetrazine based high energy density materials[)].Chi nese Journal of Energetic Materials( Hanneng Cailiao). 2017 组选择对密度的计算结果影响较小。基于力场的三维 晶体堆积方法预测的密度一般会比实验值略高,晶体结 [4] Murray JS. Brinck T. Application of calculated local surface ionization energies to chemical reactivity [] Journal of Molec 构搜索要面对巨大的状态空间和高度复杂的能量曲面 发展高效可靠的晶体结构预测方法一直是亟待解决的 5]Murray S Brinck T. Grice E et al. Correlations between mo 问题。经验公式法虽然直接但对新发现的分子或某些 lecular electrostatic potential and some experimentally-based indices of reactivity [ ]. Journal of Molecular Structure, 1992 特殊结构的适用性需验证,QSPR预测含能晶体密度的 模型精度需要提高。以下几种方案结合量子物理、计算6] Kim C K, Lee K a, Hyun K H,etal. Prediction of physico 化学、人工智能三者的优势,可显著提高含能材料密度 chemical properties of organic molecules using van der Waals surface electrostatic potentials [] Journal of Computational 计算的效率和可靠性,实现含能材料的高通量设计 Chemistry,2004,25(16):2073-2079 (1)将群体智能搜索与密度泛函理论的色散修正[7] Hohenberg P, Kohn W. In homogeneous electron gas[)]. Phys ica/ Review,1964,136(3):B864-B87 和机器学习势函数方法结合,进行含能晶体的结构和 [8] Qiu L, Xiao H, Gong X, et al. Crystal density 密度预测:先通过高斯过程回归得到势函数,用于势能 nitramines based on quantum chemistry[)]. Jour ous Materials,2007,141(1):280-288 面快速搜索,提高效率;采用全空间密网格初筛保证初[91 Gregory JO B. Joshua D H. Yonaton NH. Predicting molecu 筛样本的完备性,采用刚性分子堆积技术或晶体配位 lar crystal properties from first principles finite-temperature 搜索技术减小自由度,并适用于多种成键类型,提高成 thermochemistry to NMR crystallography [)]. Account of 功率;通过自由能表面全局最小化搜索得到备选晶体t101e8BM.Hae. Byrd E f. Accurate predictions of crystal 结构后,再用DFT-D方法对晶体结构和密度精细计 densities using quantum mechanical molecular volumes [)] 算,提高置信度 Journal of Physical Chemistry A, 2007. 111(42):10874- (2)定量构效关系与人工智能方法结合,可以方11 Politzer p,Mani. Murray Is.eta. electrostatic inter 便准确的预测含能晶体密度:建立神经网络模型描述 action correction for improved crystal density prediction [) 复杂的相关关系,结合数据挖掘,可以进一步提高模型 Molecular Physics, 2009, 107(19): 2095-2101 [12]Politzer P, Martinez J, Murray JS. et al. An electrostatic cor- 的代表性和预测精度;该方案需要高质量的含能材料 rection for improved crystal density predictions of energeti 结构与性能的实验数据做支持,为分子设计和结构-性 ionic compounds [)]. Molecular Physics, 2010, 108(10) 能关系的预测提供训练依据 l391-1396 [13]杜明燃.炸药爆轰参数理论计算方法及高能分子设计[D].安 (3)第一原理和机器学习结合(QMML方法) 徽:安徽理工大学,2016 通过高通量集成计算产生大批量的含能分子及复合体 DU Ming- ran. The research on theoretical calculations of deto- nation properties and design of high energy density compounds 系,将相关结构和性质的计算结果形成数据库,用于分 [D]. Anhui: Anhui University of Science and Technology 子表面静电势修正,基于机器学习动态更新样本结构 2016. 方面通过开发专用基函数和高精度赝势,提高第 [14] Hariji S, Uttama M, Ravi SS. Computational studies on nitro xy-1. 2. 4-triazole [l]. Journal of Ener 原理大规模并行计算能力,进行结构构造和性能预估; ic Materials,2012,30(3):265-281 另一方面,高通量计算与数据挖掘融合,可能发现隐藏15]MaP, liang I Ch, Zhu Sh C. Synthesis. XRD and DFT studies 的材料“基因”-从原子排列到相结构再到显微组织形 a novel cocrystal energetic perchlorate amine salt: methy amine triethylenediamine triperchlorate [)]. Combustion, Ex 成,最后到材料宏观性能与使用寿命之间的关系。 olosion, and Shock Waves. 2017, 53(3):319-328 Chinese Joumal of Energetic Materials. Vol28, No 1, 2020(1-12) 含能材料 www.energetic-materials.orgcn

Chinese Journal of Energetic Materials，Vol.28, No.1, 2020（1-12） 含能材料 www.energetic-materials.org.cn 王丽莉，熊鹰，谢炜宇，牛亮亮，张朝阳 5 展 望 综上所述，每种密度计算方法在拥有诸多优势的 同时都存在一定的局限性，如基于分子体积的方法在 CHNO 类含能材料中应用比较广泛，密度计算的准确 度取决于能否正确的描述分子间和分子内相互作用、 空间位阻效应、复杂的环结构以及致爆基，如 NO2、 NH2、OH、N3、ONO2、NHNO2 等 对 密 度 的 影 响 。 MSEP 优于 IED 的计算结果，一般会略微低估密度，基 组选择对密度的计算结果影响较小。基于力场的三维 晶体堆积方法预测的密度一般会比实验值略高，晶体结 构搜索要面对巨大的状态空间和高度复杂的能量曲面， 发展高效可靠的晶体结构预测方法一直是亟待解决的 问题。经验公式法虽然直接但对新发现的分子或某些 特殊结构的适用性需验证，QSPR 预测含能晶体密度的 模型精度需要提高。以下几种方案结合量子物理、计算 化学、人工智能三者的优势，可显著提高含能材料密度 计算的效率和可靠性，实现含能材料的高通量设计： （1）将群体智能搜索与密度泛函理论的色散修正 和机器学习势函数方法结合，进行含能晶体的结构和 密度预测：先通过高斯过程回归得到势函数，用于势能 面快速搜索，提高效率；采用全空间密网格初筛保证初 筛样本的完备性，采用刚性分子堆积技术或晶体配位 搜索技术减小自由度，并适用于多种成键类型，提高成 功率；通过自由能表面全局最小化搜索得到备选晶体 结构后，再用 DFT⁃D 方法对晶体结构和密度精细计 算，提高置信度。 （2）定量构效关系与人工智能方法结合，可以方 便准确的预测含能晶体密度：建立神经网络模型描述 复杂的相关关系，结合数据挖掘，可以进一步提高模型 的代表性和预测精度；该方案需要高质量的含能材料 结构与性能的实验数据做支持，为分子设计和结构⁃性 能关系的预测提供训练依据。 （3）第一原理和机器学习结合（QM\ML 方法）： 通过高通量集成计算产生大批量的含能分子及复合体 系，将相关结构和性质的计算结果形成数据库，用于分 子表面静电势修正，基于机器学习动态更新样本结构， 一方面通过开发专用基函数和高精度赝势，提高第一 原理大规模并行计算能力，进行结构构造和性能预估； 另一方面，高通量计算与数据挖掘融合，可能发现隐藏 的材料“基因”⁃从原子排列到相结构再到显微组织形 成，最后到材料宏观性能与使用寿命之间的关系。 参考文献： ［1］ Bader R F W，Carroll M T，Cheeseman J R，et al. Properties of atoms in molecule：atomic volumes［J］. Journal of American Chemical Society，1987，109（26）：7968-7979. ［2］ Rahhi S，Hari J S，Sengupta S K. Computational studies on 1，2，4⁃Triazolium⁃based salts as energetic materials［J］. Journal of Chemical Sciences，2015，127（6）：1099-1107. ［3］ 张驰，陈沫，陈潮，等 . 稠环类 1，2，4，5⁃四嗪衍生物结构和性能 的理论研究［J］. 含能材料，2017，25（4）：273-281. ZHANG Chi， CHEN Mo， CHEN Chao， et al. Theoretical study on structure and properties of polycyclic derivatives of 1，2，4，5⁃Tetrazine based high energy density materials［J］. Chi⁃ nese Journal of Energetic Materials（Hanneng Cailiao），2017， 25（4）：273-281. ［4］ Murray J S，Brinck T. Application of calculated local surface ionization energies to chemical reactivity［J］. Journal of Molec⁃ ular Structure，1992，255：271-281. ［5］ Murray J S，Brinck T，Grice E，et al. Correlations between mo⁃ lecular electrostatic potential and some experimentally⁃based indices of reactivity［J］. Journal of Molecular Structure，1992， 256：29-45. ［6］ Kim C K，Lee K A，Hyun K H，et al. Prediction of physico⁃ chemical properties of organic molecules using van der Waals surface electrostatic potentials［J］. Journal of Computational Chemistry，2004，25（16）：2073-2079. ［7］ Hohenberg P，Kohn W. In homogeneous electron gas［J］. Phys⁃ ical Review，1964，136（3）：B864-B871. ［8］ Qiu L，Xiao H，Gong X，et al. Crystal density predictions for nitramines based on quantum chemistry［J］. Journal of Hazard⁃ ous Materials，2007，141（1）：280-288. ［9］ Gregory J O B，Joshua D H，Yonaton N H. Predicting molecu⁃ lar crystal properties from first principles：finite⁃temperature thermochemistry to NMR crystallography ［J］. Account of Chemical Research，2016，49（11）：2501-2508. ［10］ Rice B M，Hare J J，Byrd E F. Accurate predictions of crystal densities using quantum mechanical molecular volumes［J］. Journal of Physical Chemistry A，2007，111（42）：10874- 10879. ［11］ Politzer P，Martinez J，Murray J S，et al. An electrostatic inter⁃ action correction for improved crystal density prediction［J］. Molecular Physics，2009，107（19）：2095-2101. ［12］ Politzer P，Martinez J，Murray J S，et al. An electrostatic cor⁃ rection for improved crystal density predictions of energetic ionic compounds［J］. Molecular Physics，2010，108（10）： 1391-1396. ［13］ 杜明燃 . 炸药爆轰参数理论计算方法及髙能分子设计［D］. 安 徽：安徽理工大学，2016. DU Ming⁃ran. The research on theoretical calculations of deto⁃ nation properties and design of high energy density compounds ［D］. Anhui：Anhui University of Science and Technology， 2016. ［14］ Hariji S，Uttama M，Ravi S S. Computational studies on nitro derivatives of 1⁃hydroxy⁃1，2，4⁃triazole［J］. Journal of Energet⁃ ic Materials，2012，30（3）：265-281. ［15］ Ma P，Jiang J Ch，Zhu Sh G. Synthesis，XRD and DFT studies of a novel cocrystal energetic perchlorate amine salt：methyl⁃ amine triethylenediamine triperchlorate［J］. Combustion，Ex⁃ plosion，and Shock Waves，2017，53（3）：319-328. 8

含能晶体密度预测的研究进展 [16]SHU Yuan-jie. LI Hua-rong. XIONG Ying, et al. Some prob lammonium 5. 5'-bistetrazole-1. 1'-diolate (TKX-50)[)] lems in theoretical design of energetic materials [)].Chinese Cryst Eng Comm,2016,18(13):2258-2267 Journal of Energetic Materials(Hanneng Cailiao), 2013, 21 [30] MENG Li-ya. LU Zhi-peng, MA Yu, et al. Enhanced interno- (2):166-172. lecular hydrogen bonds facilitating the highly dense packing of [17] Rice B M, Byrd E F. Evaluation of electrostatic descriptors for energetic hydroxylammonium Salts [). Crystal Growth &De predicting crystalline density [] Journal of Computational (12):7231-7239. Chemistry,2013.34(25):2146-2151 [31] Sorescu D C, Rice B M. Theoretical predictions of energetic [18]张君君,申程,王鹏程,等.N-1,4,6-三硝基六氢咪唑[4,5-d] molecular crystals at ambient and hydrostatic compression 咪唑-2(1H)-亚硝胺的合成与性能[门.含能材料,2017,25 conditions using dispersion corrections to conventional density (5):391-395. functionals(DFT-D)[). Journal Physical Chemistry C, 2010 ZHANG Jun-jun, SHEN Cheng, WANG Peng-cheng.et 14(14):6734-6748 Synthesis and property prediction of N(1, 4. 6-trinitrohexahy- [32] ZHNAG Lei. JIANG Seng-li, YU Yi, et al. Phase transition in droimid-azo[4. 5-d] imidazole-2(1 H)-ylidene) nitramide [) octahydro-1.3. 5. 7-tetranitro-13 5. 7-tetrazocine( HMX)un- Chinese Journal of Energetic Materials( Hanneng Cailiao) der static compression: an application of frist-principles meth 2017,25(5):391-395 od specialized for CHNO solid explosives [)]. The Journal of [19] Kim CK. Soo G C. Chang KK. et al. Prediction of densities Physical Chemistry B. 2016, 120(44):11510-11522 for solid energetic molecules with molecular surface electro-[33]余一,张蕾,姜胜利,等.TATB在高温下分解的动力学模拟[门 static potentials [I]. Journal of Computational Chemistry, 高压物理学报,2018,32(2):1-8. 2008,29(11):1818-1824 YU Yi. ZHANG Lei. JIANG Sheng-li. et al. Decomposition of [20] Nirwan A, Devi A, Ghule V D. Assessment of density predic TATB at high temperature using ab initio molecular dynamics on methods based on molecular surface electrostatic potential [)]. Chinese Journal of High Pressure Physics, 2018, 32(2) LI]. Journal of Molecular Modeling, 2018. 24(7): 166/ 166/11 [34] LIN He, CHEN Jian-fu. CUI Yu-ming, et al. A DFT-D study on [21] Ammon H L, Mitchell S. A new atom functional group volume structural. electronic, thermodynamic, and mechanical pr additivity data base for the calculation of the crystal densities erties of HMX/MPNO cocrystal under high pressure[l].Journ. of C. H. N.O and F-containing compounds []. Propellants, of Energetic Materials, 2017, 35(2):157-171 Explosives, Pyrotechnics, 1998, 23(5): 260-265 [35] Chan H CS, Kendrick ), Neumannc M A, et al. Towards ab [22]Ammon H L. Updated atom/functional group and atom_code itio screening of co-crystal formation through lattice energy volume additivity parameters for the calculation of crystal den calculations and crystal structure prediction of nicotinamide sities of single molecules, organic salts. and multi-fragment isonicotinamide. picolinamide and paracetamol multi-compo- materials containing H, C, B, N.O. F. S, P. Cl, Br. nt crystals [] Cryst Eng Comm, 2013, 15(19):3799 LI]. Propellants, Explosives, Pyrotechnics, 2008. 33(2): 92-108 [36] Panini P, Mohan T P. Gangwar U. et al. Quantitative crystal [23]Ammon H L New atom/functional group volume additivity da- tructure analysis of 1, 3, 4-thiadiazole derivatives [] Cryst ta bases for the calculation of the crystal densities of C-, H Engcomn,2013,15(22):4549-4564 N-,O-,F-,S- Cl-and Br-containing compounds [)]. [37] Podeszwa R, Rice B M. Szalewicz K. Predicting structure of Structural Chemistry.2001,12(3/4):205-212. molecular crystals from first principles[]. Physical Review Let- [24]Mathieu D. MATEO: A software package for the molecular de- ters.2008,101(11):115503/1-115503/4 ign of energetic materials[)]. Journal of Hazardous Materials, [38] Dzyabchenko A V, Pivina T S, Arnautova E A Prediction of 2010,176(1-3) tructure and density for organic nitramines[]. Journal of Me [25] Beaucamp S, Marchet N, Mathieu D, et al. Calculation of the ular Structure,1996,378(2);67-82 crystal densities of molecular salts and hydrates using additive [39] XU Xiao-juan, ZHU Wei-hua, XIAO He-ming.Hexanitro volumes for charged groups[)]. Acta Crystallographica Section hexaazaadamantane: a new potential High-Energy-Density B,2003,59(8):498-504 compound superior to hexanitrohexaazaisowurtzitane( CL-20) [26]Beaucamp S, Bernand M A. Mathieu D, et al. Ab initio sol- LI]. Journal of Energetic Materials, 2009, 27(4): 247-262 id-state heats of formation of molecular salts from ion packing [40] Srinivasan P. Asthana S N. Kumaradhas P. Ab initio crystal nd crystal modelling: application to ammonium crystals [I] tructure and charge density distribution of a highly energetic Molecular Physics, 2004. 102(3): 253-258 2. 4-dinitrobenzoic acid molecule[]. Combustion, Explosion [27]Beaucamp S, Mathieu D, Agafonov V Optimal partition ing of and Shock es,2013,49(6):712-722 olden IR, Du Z Y. Ammon H L, et al. Structure predictions crystal volumes [)]. Acta Crystallographica Section B,2007 with MOLPACK, WMIN and PMIN [C]. Universal of Maryl 63(2):277-284 USA.2009 [28] JUN Yu-fan. ZHAO Yang-zheng, YAN Su, et al. Assessment [42] Kaushik D N, Gregory JO B Prediction of organic molecular of dispersion correction methods within density functional the crystal geometries from MP2-level fragment quantum mechani- ory for energetic materials[)]. Molecular Simulation, 2017.43 cal/ molecular mechanical calculations [] The Journal of Physics,2012,137(17):174106/1-174106/11 [29] MENG Li-ya, LU Zhi-peng. WEI Xian-feng, et al. Two-sided [43] LI Jin-shan, HUANG Hui, HUANG Yi-gang et al. A Theoreti- effects of strong hydrogen bonding on the stability of dihydro- cal study of highly nitrated azacubanes[]. Journal of Energetic CHINESE JOURNAL OF ENERGETIC MATERIALS 含能材料 2020年第28卷第1期(1-12

CHINESE JOURNAL OF ENERGETIC MATERIALS 含能材料 2020 年 第 28 卷 第 1 期 （1-12） 含 能 晶 体 密 度 预 测 的 研 究 进 展 ［16］ SHU Yuan⁃jie，LI Hua⁃rong，XIONG Ying，et al. Some prob⁃ lems in theoretical design of energetic materials［J］. Chinese Journal of Energetic Materials（Hanneng Cailiao），2013，21 （2）：166-172. ［17］ Rice B M，Byrd E F. Evaluation of electrostatic descriptors for predicting crystalline density［J］. Journal of Computational Chemistry，2013，34（25）：2146-2151. ［18］ 张君君，申程，王鹏程，等 . N⁃1，4，6⁃三硝基六氢咪唑［4，5⁃d］ 咪 唑⁃2（1H）⁃亚 硝 胺 的 合 成 与 性 能［J］. 含 能 材 料 ，2017，25 （5）：391-395. ZHANG Jun⁃jun，SHEN Cheng，WANG Peng⁃cheng，et al. Synthesis and property prediction of N（⁃ 1，4，6⁃trinitrohexahy⁃ droimid⁃azo［4，5⁃d］imidazole⁃2（1H）⁃ylidene） nitramide［J］. Chinese Journal of Energetic Materials（Hanneng Cailiao）， 2017，25（5）：391-395. ［19］ Kim C K，Soo G C，Chang K K，et al. Prediction of densities for solid energetic molecules with molecular surface electro⁃ static potentials ［J］. Journal of Computational Chemistry， 2008，29（11）：1818-1824. ［20］ Nirwan A，Devi A，Ghule V D. Assessment of density predic⁃ tion methods based on molecular surface electrostatic potential ［J］. Journal of Molecular Modeling，2018，24（7）：166/1- 166/11. ［21］ Ammon H L，Mitchell S. A new atom functional group volume additivity data base for the calculation of the crystal densities of C，H，N，O and F⁃containing compounds［J］. Propellants， Explosives，Pyrotechnics，1998，23（5）：260-265. ［22］ Ammon H L. Updated atom/functional group and atom_code volume additivity parameters for the calculation of crystal den⁃ sities of single molecules，organic salts，and multi⁃fragment materials containing H，C，B，N，O，F，S，P，Cl，Br，and I ［J］. Propellants， Explosives， Pyrotechnics， 2008， 33（2）： 92-108. ［23］ Ammon H L. New atom/functional group volume additivity da⁃ ta bases for the calculation of the crystal densities of C⁃，H⁃， N⁃，O⁃，F⁃，S⁃，P⁃，Cl⁃，and Br⁃containing compounds［J］. Structural Chemistry，2001，12（3/4）：205-212. ［24］ Mathieu D. MATEO：A software package for the molecular de⁃ sign of energetic materials［J］. Journal of Hazardous Materials， 2010，176（1-3）：313-322. ［25］ Beaucamp S，Marchet N，Mathieu D，et al. Calculation of the crystal densities of molecular salts and hydrates using additive volumes for charged groups［J］. Acta Crystallographica Section B，2003，59（8）：498-504. ［26］ Beaucamp S，Bernand M A，Mathieu D，et al. Ab initio sol⁃ id⁃state heats of formation of molecular salts from ion packing and crystal modelling：application to ammonium crystals［J］. Molecular Physics，2004，102（3）：253-258. ［27］ Beaucamp S，Mathieu D，Agafonov V. Optimal partitioning of molecular properties into additive contributions：the case of crystal volumes［J］. Acta Crystallographica Section B，2007， 63（2）：277-284. ［28］ JUN Yu⁃fan，ZHAO Yang⁃zheng，YAN Su，et al. Assessment of dispersion correction methods within density functional the⁃ ory for energetic materials［J］. Molecular Simulation，2017，43 （7）：568-574. ［29］ MENG Li⁃ya，LU Zhi⁃peng，WEI Xian⁃feng，et al. Two⁃sided effects of strong hydrogen bonding on the stability of dihydrox⁃ ylammonium 5，5′⁃bistetrazole⁃1，1′⁃diolate （TKX⁃50）［J］. Cryst Eng Comm，2016，18（13）：2258-2267. ［30］ MENG Li⁃ya，LU Zhi⁃peng，MA Yu，et al. Enhanced intermo⁃ lecular hydrogen bonds facilitating the highly dense packing of energetic hydroxylammonium Salts［J］. Crystal Growth & De⁃ sign，2016，16（12）：7231-7239. ［31］ 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］. Journal Physical Chemistry C，2010， 114（14）：6734-6748. ［32］ ZHNAG Lei，JIANG Seng⁃li，YU Yi，et al. Phase transition in octahydro⁃1，3，5，7⁃tetranitro⁃1，3，5，7⁃tetrazocine（HMX）un⁃ der static compression：an application of frist⁃principles meth⁃ od specialized for CHNO solid explosives［J］. The Journal of Physical Chemistry B，2016，120（44）：11510-11522. ［33］ 余一，张蕾，姜胜利，等 . TATB 在高温下分解的动力学模拟［J］. 高压物理学报，2018，32（2）：1-8. YU Yi，ZHANG Lei，JIANG Sheng⁃li，et al. Decomposition of TATB at high temperature using ab initio molecular dynamics ［J］. Chinese Journal of High Pressure Physics，2018，32（2）： 1-8. ［34］ LIN He，CHEN Jian⁃fu，CUI Yu⁃ming，et al. A DFT⁃D study on structural，electronic，thermodynamic，and mechanical prop⁃ erties of HMX/MPNO cocrystal under high pressure［J］. Journal of Energetic Materials，2017，35（2）：157-171. ［35］ Chan H C S，Kendrick J，Neumannc M A，et al. Towards ab initio screening of co⁃crystal formation through lattice energy calculations and crystal structure prediction of nicotinamide， isonicotinamide，picolinamide and paracetamol multi⁃compo⁃ nent crystals［J］. Cryst Eng Comm，2013，15（19）：3799- 3807. ［36］ Panini P，Mohan T P，Gangwar U，et al. Quantitative crystal structure analysis of 1，3，4⁃thiadiazole derivatives［J］. Cryst EngComm，2013，15（22）：4549-4564. ［37］ Podeszwa R，Rice B M，Szalewicz K. Predicting structure of molecular crystals from first principles［J］. Physical Review Let⁃ ters，2008，101（11）：115503/1-115503/4. ［38］ Dzyabchenko A V，Pivina T S，Arnautova E A. Prediction of structure and density for organic nitramines［J］. Journal of Mo⁃ lecular Structure，1996，378（2）：67-82. ［39］ XU Xiao⁃juan， ZHU Wei⁃hua， XIAO He⁃ming. Hexanitro⁃ hexaazaadamantane： a new potential High⁃Energy⁃Density compound superior to hexanitrohexaazaisowurtzitane（CL⁃20） ［J］. Journal of Energetic Materials，2009，27（4）：247-262. ［40］ Srinivasan P，Asthana S N，Kumaradhas P. Ab initio crystal structure and charge density distribution of a highly energetic 2，4⁃dinitrobenzoic acid molecule［J］. Combustion，Explosion， and Shock Waves，2013，49（6）：712-722. ［41］ Holden J R，Du Z Y，Ammon H L，et al. Structure predictions with MOLPACK，WMIN and PMIN［C］. Universal of Maryl， USA，2009. ［42］ Kaushik D N，Gregory J O B. Prediction of organic molecular crystal geometries from MP2⁃level fragment quantum mechani⁃ cal/molecular mechanical calculations ［J］. The Journal of Chemical Physics，2012，137（17）：174106/1-174106/11. ［43］ LI Jin⁃shan，HUANG Hui，HUANG Yi⁃gang，et al. A Theoreti⁃ cal study of highly nitrated azacubanes［J］. Journal of Energetic 9

王丽莉,熊鹰,谢炜宇,牛亮亮,张朝阳 Materials,2008,26(4):230-245. [60]Simoes P N. Pedroso L M. Matos B A M, et al. Crystal and [44] Ghule V D. Radhakrishnan S, Jadhav P M, et al. Quantum molecular structure of 4. 6-Bis (nitroimino)-1. 3. 5-triazi- chemical studies on energetic Azo-Bridged azoles [)]. Journal an-2-one: theoretical and X-ray studies [] Journal of Physical of Energetic Materials, 2013, 31(1): 35-48. Chemistry A.2007,11101):150-158 [45] Magnus B, Carl C. A validation study of the general [61] Pickard C I. Needs R I. High-pressure phases of silane [ force field applied to energetic molecular crystals [l] Physical Review Letters,2006,97(4):045504/1-045504/4 of Energetic Materials. 2016, 34(1): 62-75. [62]Zilka M. Dudenko D V, Hughes C E. et al. Ab initio random [46] Ravi P, Gore G M. Tewari Surya P. et al. Theoretical studies structure searching of organic molecular solids: assessment on aminoand Methyl-Substituted trinitrodiazoles [)]. Journal of and validation against experimental data [] Physical chemis- Energetic Materials, 2011. 29(3): 209-227. ry chemical Physics,2017,19(38):25949-25960 [47] Han W, Habgood M, Parker), et al. Crystal structure deter- [63] Lin T J, Hsing CR, Wei M, et al. Structure prediction of the mination by combined synchrotron powder X-ray diffraction lid forms of methanol: an ab initio random structure search and crystal structure prediction: 1: 1 L-ephedrine D-tartrate ing approach [)]. Physical chemistry chemical Physics,2016 [刀. Cryst Eng Comm,2013.15(10):1853-1859 18(4):2736-2746 [48] Ghule V D. Nirwan A Role of forcefield in density prediction [64] Oganov A R, Lyakhow AO. Valle M. How evolutionary crys- for CHNO explosives [)]. Structural Chemistry, 2018. 29(5) al structure prediction works, and why[). Account of Chemi- cal research,2011.44(3):227-237 [49]王喜军,赵立峰,孙淮.TNAz的第一原理分子力学力场[]含[65] Qiang Z, Artem R C, Andriy O L. Evolutionary metadynam 能材料,200715(6):616-621 ics: a novel method to predict crystal structures[)]. CrystEng WANG Xi-jun, ZHAO Li-feng. SUN Huai. The frist principle Comm,2012,14(10):3596-3601 orce field for energetic material TNAZ[门]. Chinese Journal of[66]彭军辉,曾庆丰,谢聪伟,等.HC体系的高压结构预测及电子 Energetic Materials(Hanneng Cailiao). 2007. 15(6) 性质第一性原理模拟[冂].物理学报,2015,64(23):236102/1 616-621 236102/10. [50] XUE Xiang-gui. WEN Yu-shi, ZHANG Chao-yang. Early de PENG Jun-hui, ZENG Qing-feng. XIE Cong-wei, et al. cay mechanism of shocked E-CL-20: a molecular dynamics High-pressure structure prediction of Hf-C system and simulation study [] Joumal Physical Chemistry C.2016, 120 first-principle simulation of their electronic properties [)].Acta (38):21169-21177 Physics Sinica,2015,64(23):236102/1-236102/10 [51]Holden JR, Du Z Y. Ammon H L Prediction of possible crys- [67] Arapan S Nieves P, Cuesta-L6pez S A high-throughput explo- tal structures for C-, H-.N-.O-, and F-containing organ- ration of magnetic materials by using structure predicting meth- ic compounds[)]. Journal of Computational Chemistry, 1993 ods[)]. Journal of Applied Physics, 2018, 123(8):083904/ 14(4):422-437 [52] Reilly A M. Cooper R I, Adjiman CS. et al. Report on the [68] Zhang M G. Wang H, Wang H B, et al. Frist-Principles pre sixth blind test of organic crystal structure prediction methods diction on the high-pressure structures of transition metal dibo [1]. Acta Crystallographica Section B,2016,72(4):439-459. des( TMB. TM=Sc, Ti, Y, Zr)[]. Inorganic chemistry, 2010 [53]Kirkpatrick S, Gelatt C D, Vecchi M P, et al. Optimization by 49(15):6859-6866 simulated annealing[ ]].Science, 1983. 220: 671-680 [69]王洪波.高压下含氮高能量高密度材料的设计与合成[D].吉 [54]Metropolis N, Rosenbluth A. Rosenbluth M, et al. Equation of 林:吉林大学,2012 state calculations by fast computing mathies []. Journal of WANG Hong-bo. Theoretical design and experimental synth Chemical Physics,1953,56(21):1087-109 sis of nitrogen containing high energy density material at high [55] Yin 1, Chaitanya K, Ju X H Bromination and cyanation for im pressure[ D]. Jilin: Jilin University, 2012 proving electron transport performance of anthra- tetrathio-[70]吕健 CALYPSO晶体和团簇结构预测方法及应用[D].吉林:吉 phene []. Journal of Materials Research, 2016, 31(3) 林大学,2013. 337-347 LV Jian. CALYPSO crystal and cluster structure prediction [56]ZHAO Guo-zheng, LU Ming Symmetry. conjugated system method and application[D]. Jilin: Jilin University, 2013 and N-oxides in 2. 4. 6-trinitro-1. 3. 5-triazine-1.3. 5-trioxide [71] Wang Y C, Lv J. Zhu L, et al. Materials discovery via CALYP a new design concept for energetic oxidizers[]. Journal of En- SO methodology [)]. Joumal of physics: condensed matter. ergetic Materials, 2017, 35(1):63-76 2015,27(20):203203-203217 [57]刘英晢,来蔚鹏,尉涛,等.全氮材料基础性能理论研究:I.晶[72]王晓丽,李建福,陈丽.基于 CALYPSO方法的新型高能量密度 体密度预测[].含能材料,2017,25(2):100-105 材料设计[门].科学通报,2015,60(27):2608-2615 LIU Ying-zhe, LAI Wei-peng, YU Tao, et al. Theoretical inves- WANG Xiao-li, LI Jian-fu, CHEN Li. The design of new high tigations on fundamental properties of All-Nitrogen materials energy density materials based on CALYPSO method [].Chin 1. prediction of crystal densities[)]. Chinese journal of Energet- Sci Bull,2015,60(27):2608-2615 ic Materials( Hanneng Cailiao).2017. 25(2): 100-105 [73]王彦超,吕健,马琰铭. CALYPSO结构预测方法[J].科学通 [58] Moxnes )F. Hansen F K, Jensen T L, et al. A computational 报,2015.60(27):2580-2587 study of density of some high energy molecules [)]. Prop WANG Yan-chao, LV J. MA Yan-ming. CALYPSO structure od [j [59] Meyer R, Kohler 1. Homburg A. Explosives, 6 th Edition [M] Germany: Wiley-VCH, Weinheim. 2007. 11 [74]宗朝霞,汤宏胜,贺曼,等.基于遗传算法的支持向量机预测含 Chinese Joumal of Energetic Materials. Vol28, No 1, 2020(1-12) 含能材料 www.energetic-materials.orgcn

Chinese Journal of Energetic Materials，Vol.28, No.1, 2020（1-12） 含能材料 www.energetic-materials.org.cn 王丽莉，熊鹰，谢炜宇，牛亮亮，张朝阳 Materials，2008，26（4）：230-245. ［44］ Ghule V D，Radhakrishnan S，Jadhav P M，et al. Quantum chemical studies on energetic Azo⁃Bridged azoles［J］. Journal of Energetic Materials，2013，31（1）：35-48. ［45］ Magnus B，Carl C. A validation study of the general Amber force field applied to energetic molecular crystals［J］. Journal of Energetic Materials，2016，34（1）：62-75. ［46］ Ravi P，Gore G M，Tewari Surya P，et al. Theoretical studies on aminoand Methyl⁃Substituted trinitrodiazoles［J］. Journal of Energetic Materials，2011，29（3）：209-227. ［47］ Han W，Habgood M，Parker J E，et al. Crystal structure deter⁃ mination by combined synchrotron powder X⁃ray diffraction and crystal structure prediction：1：1 L⁃ephedrine D⁃tartrate3 ［J］. Cryst Eng Comm，2013，15（10）：1853-1859. ［48］ Ghule V D，Nirwan A. Role of forcefield in density prediction for CHNO explosives［J］. Structural Chemistry，2018，29（5）： 1-8. ［49］ 王喜军，赵立峰，孙淮 . TNAZ 的第一原理分子力学力场［J］. 含 能材料，2007，15（6）：616-621. WANG Xi⁃jun，ZHAO Li⁃feng，SUN Huai. The frist principle force field for energetic material TNAZ［J］. Chinese Journal of Energetic Materials （Hanneng Cailiao） ， 2007， 15 （6） ： 616-621. ［50］ XUE Xiang⁃gui，WEN Yu⁃shi，ZHANG Chao⁃yang. Early de⁃ cay mechanism of shocked ε⁃CL⁃20：a molecular dynamics simulation study［J］. Journal Physical Chemistry C. 2016，120 （38）：21169-21177. ［51］ Holden J R，Du Z Y，Ammon H L. Prediction of possible crys⁃ tal structures for C─，H─，N─，O─，and F⁃containing organ⁃ ic compounds［J］. Journal of Computational Chemistry，1993， 14（4）：422-437. ［52］ Reilly A M，Cooper R I，Adjiman C S，et al. Report on the sixth blind test of organic crystal structure prediction methods ［J］. Acta Crystallographica Section B，2016，72（4）：439-459. ［53］ Kirkpatrick S，Gelatt C D，Vecchi M P，et al. Optimization by simulated annealing［J］. Science，1983，220：671-680. ［54］ Metropolis N，Rosenbluth A，Rosenbluth M，et al. Equation of state calculations by fast computing mathines［J］. Journal of Chemical Physics，1953，56（21）：1087-1092. ［55］ Yin J，Chaitanya K，Ju X H. Bromination and cyanation for im ⁃ proving electron transport performance of anthra⁃tetrathio⁃ phene［J］. Journal of Materials Research， 2016， 31（3）： 337-347. ［56］ ZHAO Guo⁃zheng，LU Ming. Symmetry，conjugated system， and N⁃oxides in 2，4，6⁃trinitro⁃1，3，5⁃triazine⁃1，3，5⁃trioxides： a new design concept for energetic oxidizers［J］. Journal of En⁃ ergetic Materials，2017，35（1）：63-76. ［57］ 刘英哲，来蔚鹏，尉涛，等 . 全氮材料基础性能理论研究：Ⅰ. 晶 体密度预测［J］. 含能材料，2017，25（2）：100-105. LlU Ying⁃zhe，LAI Wei⁃peng，YU Tao，et al. Theoretical inves⁃ tigations on fundamental properties of All⁃Nitrogen materials： 1. prediction of crystal densities［J］. Chinese journal of Energet⁃ ic Materials（Hanneng Cailiao），2017，25（2）：100-105. ［58］ Moxnes J F，Hansen F K，Jensen T L，et al. A computational study of density of some high energy molecules［J］. Propel⁃ lants，Explosives，Pyrotechnics，2017，42（2）：204-212. ［59］ Meyer R，Köhler J，Homburg A. Explosives，6th Edition［M］. Germany：Wiley⁃VCH，Weinheim，2007，11. ［60］ Simöes P N，Pedroso L M，Matos B A M，et al. Crystal and molecular structure of 4，6⁃Bis（nitroimino）⁃1，3，5⁃triazi⁃ nan⁃2⁃one：theoretical and X⁃ray studies［J］. Journal of Physical Chemistry A，2007，111（1）：150-158. ［61］ Pickard C J，Needs R J. High⁃pressure phases of silane［J］. Physical Review Letters，2006，97（4）：045504/1-045504/4. ［62］ Zilka M，Dudenko D V，Hughes C E，et al. Ab initio random structure searching of organic molecular solids：assessment and validation against experimental data［J］. Physical chemis⁃ try chemical Physics，2017，19（38）：25949-25960. ［63］ Lin T J，Hsing C R，Wei C M，et al. Structure prediction of the solid forms of methanol：an ab initio random structure search⁃ ing approach［J］. Physical chemistry chemical Physics，2016， 18（4）：2736-2746. ［64］ Oganov A R，Lyakhow A O，Valle M. How evolutionary crys⁃ tal structure prediction works，and why［J］. Account of Chemi⁃ cal Research，2011，44（3）：227-237. ［65］ Qiang Z，Artem R O，Andriy O L. Evolutionary metadynam⁃ ics：a novel method to predict crystal structures［J］. CrystEng⁃ Comm，2012，14（10）：3596-3601. ［66］ 彭军辉，曾庆丰，谢聪伟，等 . Hf⁃C 体系的高压结构预测及电子 性质第一性原理模拟［J］. 物理学报，2015，64（23）：236102/1- 236102/10. PENG Jun⁃hui， ZENG Qing⁃feng， XIE Cong⁃wei， et al. High⁃pressure structure prediction of Hf⁃C system and first⁃principle simulation of their electronic properties［J］. Acta Physics Sinica，2015，64（23）：236102/1-236102/10. ［67］ Arapan S，Nieves P，Cuesta⁃López S. A high⁃throughput explo⁃ ration of magnetic materials by using structure predicting meth⁃ ods［J］. Journal of Applied Physics，2018，123（8）：083904/1- 083904/8. ［68］ Zhang M G，Wang H，Wang H B，et al. Frist⁃Principles pre⁃ diction on the high⁃pressure structures of transition metal dibo⁃ rides（TMB2，TM=Sc，Ti，Y，Zr）［J］. Inorganic chemistry，2010， 49（15）：6859-6866. ［69］ 王洪波 . 高压下含氮高能量高密度材料的设计与合成［D］. 吉 林：吉林大学，2012. WANG Hong⁃bo. Theoretical design and experimental synthe⁃ sis of nitrogen containing high energy density material at high pressure［D］. Jilin：Jilin University，2012. ［70］ 吕健 .CALYPSO 晶体和团簇结构预测方法及应用［D］. 吉林：吉 林大学，2013. LV Jian. CALYPSO crystal and cluster structure prediction： method and application［D］. Jilin：Jilin University，2013. ［71］ Wang Y C，Lv J，Zhu L，et al. Materials discovery via CALYP⁃ SO methodology［J］. Journal of physics：condensed matter， 2015，27（20）：203203-203217. ［72］ 王晓丽，李建福，陈丽 . 基于 CALYPSO 方法的新型高能量密度 材料设计［J］. 科学通报，2015，60（27）：2608-2615. WANG Xiao⁃li，LI Jian⁃fu，CHEN Li. The design of new high energy density materials based on CALYPSO method［J］. Chin Sci Bull，2015，60（27）：2608-2615. ［73］ 王 彦 超 ，吕 健 ，马 琰 铭 . CALYPSO 结 构 预 测 方 法［J］. 科 学 通 报，2015，60（27）：2580-2587. WANG Yan⁃chao，LV J，MA Yan⁃ming. CALYPSO structure prediction method［J］. Chin Sci Bull，2015，60（27）：2580- 2587. ［74］ 宗朝霞，汤宏胜，贺曼，等 . 基于遗传算法的支持向量机预测含 能 材 料 密 度 的 研 究［J］. 计 算 机 与 应 用 化 学 ，2009，26（12）： 1529-1533. ZONG Zhao⁃xia，TANG Hong⁃sheng，HE Man，et al. Re⁃ search on QSPR for energetic materials based on genetic algo⁃ rithm support vector machine［J］. Computers and applied Chemistry，2009，26（12）：1529-1533. ［75］ Cho S G，Goh E M，Kim J K. Holographic QSAR models for es⁃ timating densities of energetic materials［J］. Bull Korean Chem Soc，2001，22（7）：775-778. ［76］ Lee A，Kim D，Kim K H，et al. Elucidation of specific aspects of dielectric constants of conjugated organic compounds：a QSPR approach［J］. Journal of Molecular Modeling，2012，18 （1）：251-256. ［77］ 高森 . 混合炸药感度理论预测研究［D］. 山西：中北大学，2011. GAO Sen. Research in sensitivity prediction of composite ex⁃ plosives［D］. Shanxi：North University of China，2011. ［78］ Xu J，Zhu L，Fang D，et al. QSPR studies of impact sensitivity of nitro energetic compounds using Three⁃Dimensional de⁃ scriptors［J］. Journal of Molecular Graphics and Modelling， 2012，36（6）：10-19. ［79］ 来蔚鹏，廉鹏，王伯周，等 . 用定量结构性质关系（QSPR）预测 芳香系炸药的密度［J］. 含能材料，2007，15（6）：626-655. LAI Wei⁃peng，LIAN Peng，WANG Be⁃zhou，et al. Prediction of density of aromatic explosives by quantitative struc⁃ ture⁃property relationship（QSPR）method［J］. Chinese Journal of Energetic Materials（Hanneng Cailiao），2007， 15（6）： 626-655. ［80］ Hwang S N，Lee J H，Lee S K，et al. Rapid and robust QSPR model for prediction of density for high energetic materials （HEMs）［C］//3th Symposia on High Energetic materials. ［81］ Li Jun⁃ling，Lu Fang⁃yun，Qin Jin⁃gui，et al. Effects of tempera⁃ ture and strain rate on the dynamic responses of three poly⁃ mer⁃bonded explosives［J］. Journal of Strain Analysis for Engi⁃ neering Design，2012，47（2）：104-112. ［82］ Mehmet S O. Artificial neural network approach to predict the electrical conductivity and density of Ag⁃Ni binary alloys［J］. Journal of Materials Processing Technology，2008，208（1-3）： 470-476. ［83］ 黄俊，周申范 . 人工神经网络法预测炸药爆速的研究［J］. 火炸药 学报，2000，23（1）：34-37. HUANG Jun，ZHOU Shen⁃fan. Predicting the detonating ve⁃ locity of explosives by artificial neural network［J］. Chinese Journal of Explosive & Propellants，2000，23（1）：34-37. ［84］ Wang D，Yuan Y，Duan S，et al. QSPR study on melting point of carbocyclic nitroaromatic compounds by multiple linear re⁃ gression and artificial neural network［J］. Chemometrics and In⁃ telligent Laboratory Systems，2015，143（1）：7-15. ［85］ Tatyana K. Detonation properties of high explosives calculated by revised Kihara⁃Hikita equation of state［C］//The Eighth Sym⁃ posium（International）on Detonation，1985：548-557. ［86］ 荔建锋，王蕾娜，陈念贻 . 人工神经网络用于 PBX 炸药装药密度 的研究［J］，火炸药学报，1997，（2）：26-29. LI Jian⁃feng，WANG Lei⁃na，CHEN Nian⁃yi. Research the den⁃ sity of PBX using aritificial neural network［J］. Chinese Journal of Explosive & Propellants，1997，20（2）：26-29. ［87］ 王国栋 . 含能材料主要性能参数的预测研究［D］. 山西：中北大 学，2007. WANG Guo⁃dong. Study on predicting properties of energetic material［D］. Shanxi：North University of China，2007. ［88］ 王国栋，刘存玉 . 神经网络在炸药晶体密度预测中的应用［J］. 火 炸药学报 .2007，30（1）：57-63. WANG Guo⁃dong，LIU Cun⁃yu. Application of artificial neural network in predicting the density of explosives［J］. Chinese Journal of Explosives & Propellants，2007，30（1）：57-63. ［89］ 章婷曦，黄俊，周申范 . 人工神经网络法预测炸药组分的色谱保 留值参数［J］. 色谱，2001，19（4）：319-322. ZHANG Ting⁃xi，HUANG Jun，ZHOU Shen⁃fan. Prediction of retention parameters of explosives by artificial neural network ［J］. Chinese Journal of Chromatography， 2001， 19（4）： 319-322. ［90］ 程娥，陈莉，葛忠学，等 . 一种改进的遗传算法在含能材料设计 中 的 研 究 与 应 用［J］. 计 算 机 与 应 用 化 学 ，2008，25（2）： 164-168. CHENG E，CHEN Li，GE Zhong⁃xue，et al. Application and re⁃ search of advanced GA in QSAR studies［J］. Computers and ap⁃ plied Chemistry，2008，25（2）：164-168. ［91］ Keshavarz M H. Prediction of densities of acyclic and cyclic ni⁃ tramines，nitrate esters and nitroaliphatic compounds for eval⁃ uation of their detonation performance［J］. Journal of Hazard⁃ ous Materials，2007，143（1-2）：437-442. ［92］ Keshavarz M H. New method for calculating densities of ni⁃ troaromatic explosive compounds［J］. Journal of Hazardous Materials，2007，145（1-2）：263-269. ［93］ Keshavarz M H，Pouretedal H R. A reliable simple method to estimate density of nitroaliphatics，nitrate esters and nitramines ［J］. Journal of Hazardous Materials， 2009， 169（1-3） ： 158-169. ［94］ Keshavarz M H，Zohari N，Seyedsadjadi S A. Validation of im⁃ proved simple method for prediction of activation energy of the thermal decomposition of energetic compounds［J］. Journal of Thermal Analysis and Calorimetry， 2013， 114 （2） ： 497-510. ［95］ Keshavarz M H，Motamedoshariati H，Moghayadnia R，et al. A new computer code for assessment of energetic materials with crystal censity，condense phase enthalpy of formation， and activation energy of thermolysis［J］. Propellants，Explo⁃ sive，Pyrotechnics，2013，38（1）：95-102. ［96］ Keshavarz M H，Rahimi R，Akbarzadeh A R. Two novel corre⁃ lations for assessment of crystal density of hazardous ionic mo⁃ lecular energetic materials using their molecular structures［J］. Fluid Phase Equilibria，2015，402（1）：1-8. ［97］ Keshavarz M H. Novel method for predicting densities of polynitro arene and polynitro heteroarene explosives in order to evaluate their detonation performance［J］. Journal of Hazard⁃ ous Materials，2009，165（2-3）：579-588. ［98］ Frem D. Theoretical studies on energetic properties of s⁃Tri⁃ azine substituted aminofurazan and aminofurazan and amino⁃ furoxan derivatives⁃high performance energetic material sys⁃ tems［J］. Combustion，Explosion，and Shock Waves，2014，50 （4）：441-446. ［99］ Keshavarz M H， Soury H，Motamedoshariati H，et al. Im⁃ proved method for prediction of density of energetic com ⁃ pounds using their molecular structure［J］. Structure Chemis⁃ try，2015，26（2）：455-466. ［100］Keshavarz M H，Jafari M，Motamedoshariati H，et al. Energet⁃ ic materials designing bench（EMDB），Version 1.0［J］. Propel⁃ 10