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

含能晶体密度预测的研究进展 文章编号: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