高通量计算与深度学习相结合的稠环含能化合物设计

王润文，杨春明，刘 建

（1. 西南科技大学计算机科学与技术学院，四川 绵阳 621000；2. 中国工程物理研究院化工材料研究所，四川 绵阳 621999）

0 引言

含能分子是一类具有高能量属性的特殊分子，广泛用于国防、航空航天和民用领域［1-4］。含能分子在研制的各环节中必须兼顾多项性能指标，包括能量［5-7］、吸湿性［8-10］、热稳定性［11-12］、化学键能［13-14］、感度［15-16］等。因此，新含能分子从最初的设计到最终的应用需要经历复杂的流程，发现真正有应用前景的新含能分子极具挑战性。20 世纪80 年代开始，计算机辅助材料设计（Compute-aided material design，CAMD）［17］开始兴起，并逐步被引入到含能材料的研究中，使含能材料的性能可以通过计算机手段来获得［18-21］。CAMD 可以通过预测含能分子的性质提前评估其是否具备合成价值，从而筛选出有价值的分子，减少了实验中的“试错”行为。在CAMD 的协助下，含能分子的研制模式已由传统的“试错法”升级为计算与实验协同的模式，即先由计算筛选出满足性质的目标含能分子，然后探索合成路线。得益于此，化学家可以避免在无价值的分子上耗费时间和资源，使含能分子的研制效率在一定程度上得到提升。

在实际应用中，含能分子设计的成功率取决于2个主要因素，即分子结构的构建方式和性能预测方法。在分子结构的构建方面，由于不同的元素组成和原子间连接存在多种可能性，CAMD 的研究对象往往是由一定数量的分子结构组成集合，也即虚拟筛选空间［22］。在含能分子的性能预测方面，一般采用量子化学计算或定量构效关系（Quantitative structure-property relationship，QSPR）模型［19，23-24］来实现。含能分子设计的本质就是基于性能预测的结果对虚拟筛选空间中的分子结构进行排序，排序靠前的分子结构将作为候选含能分子推荐给合成化学家。因此，为保证含能分子设计的“成功率”，高效的虚拟筛选空间构建方式和准确的性能预测方法是基本前提。近年来，基于机器学习的QSPR 和高通量筛选技术开始用于CAMD，并在含能分子设计中获得了广泛的应用，使含能分子的设计水平大幅提升。

Pinheiro 等［25］通过机器学习预测QM9 数据 集中分子的9 种性质，最终模型精度最高HUMO-LUMO预测精度（R2）达到了0.97。Kang 等［26］通过机器学习预测含能分子的生成焓等性质，其训练的模型的R2达到0.93。Elton 等［27］以键总数为描述符，采用Kernel Ridge 回归，得到晶体密度性质的预测模型，R2为0.74。Casey 等［28］采用分子三维电子结构为描述符，基于Deep Neural Network（DNN）得到含能材料7 种主要性质的预测模型，7 种模型的R2值均达到0.95 以上。徐雅斌等［29］以库仑矩阵为输入，基于Convolutional Neural Network（CNN）和Bi-directional Short-term Memory Network（Bi-LSTM）的融合模型，实现了对含能材料生成焓的预测，其训练出的模型的R2为0.97。Yang 等［30］基于图神经网络训练了一个密度预测模型，可以通过分子结构直接预测出晶体密度，模型精度R2值达到了0.949。Jiang 等［31］通过神经网络获得了分子的晶体的共晶可行性判据模型，在此基础上成功制备出一种新型含能共晶。Huang 等［32］采用晶体描述符和梯度提升机（XGBoost）等5 种算法获得了面向爆热、爆速、爆压、热分解温度和晶格能的机器学习模型，建立了综合评估爆轰性能和稳定性的高通量解决方案。Liu 等［33］运用高通量筛选技术设计并筛选了大量含 有CHONF 的 含 能 分 子。Song 等［34］基 于 机 器 学 习建立了分子性质预测和晶体堆积模式评估模型，从25 112 个分子构成的虚拟筛选空间中筛选出综合性优异的新型含能分子，并实现了实验室合成与测试。

上述研究工作将机器学习和高通量筛选技术用于含能分子的设计，为含能分子的高效筛选提供了解决方案，但在虚拟筛选空间构建和性能预测模型方面仍有改进的空间。对于虚拟筛选空间构建，当前的CAMD 需要研究者依据知识和经验选定基本结构单元和官能团种类后再逐一枚举分子结构。该方法以人为判断为前提，易造成虚拟筛选空间构建效率低和分子结构类型覆盖面窄等问题。对于性能预测方法，有多种方法可以实现对爆轰性能准确的预测，且与爆轰性能密切相关的密度和生成焓预测结果一般较可靠。存在的问题是部分模型基于的训练样本少、适用范围窄，尤其是化学稳定性相关的性质（比如感度）一直以来都难以被准确预测。

本研究拟将含能分子研究中广受关注的稠环类结构为研究对象，开展高通量计算与深度学习相结合的含能分子设计，并针对上述问题提出相应的解决方案。在虚拟筛选空间构建中，基于分子骨架等效密度的预筛选获得了高密度分子骨架，然后再通过分子片段组装构建了具有潜在高密度特征的虚拟筛选空间。在性能预测中，根据不同预测方法的适用性，（i）采用量子化学计算预测生成焓、X-NO2键解离能（Bond Dissociation Energy，BDE），（ii）采用深度学习模型预测晶体密度；（iii）采用爆轰性能经验公式预测爆轰性能。针对含能分子化学稳定性难以准确预测的问题，采用与化学稳定性密切相关的X-NO2BDE［56］间接评估化学稳定性。在高密度骨架预筛选的基础上构建了具有潜在高密度的虚拟筛选空间，并将高通量计算和机器学习相结合预测了密度、生成焓、爆速、爆压和X-NO2BDE，从中筛选出性能优异的候选含能分子共计6 种。

1 原理与方法

1.1 基本原理

基于CAMD 的含能分子筛选，首先需要构建由大量待筛选分子构成的虚拟筛选空间，然后采用量子化学计算或是QSPR 模型进行筛选，最终获得符合性能预期的含能分子。高效的分子筛选需要解决2 个问题：一是如何构建具有高性能潜质筛选对象，以提高分子筛选的成功率；二是如何获得精度高、泛化性好的模型以实现对性能的准确预测。本研究以近年来较热门的稠环类分子作为研究对象，基于高通量计算与深度学习相结合的方法，力图发现高密度、高爆轰性能和高化学稳定性的含能分子候选物。具体的技术原理如图1 所示。

图1 高通量计算与深度学习相结合的稠环含能分子设计技术原理Fig.1 Principle for designing fused-ring energetic compounds using high-throughput computing combined with deep learning

如图1 所示，研究首先预先筛选高密度分子骨架，然后经分子片段组装获得由潜在的高密度含能分子构成的虚拟筛选空间。在此基础上基于深度学习、量子化学计算和爆轰产物状态方程等方法预测晶体密度、生成焓、爆轰性能、BDE 和撞击感度的数值，从而由性能排序完成含能分子的筛选。虚拟筛选空间的构建由分子高密度骨架筛选和分子片段对接两步构成。首先，从收集到的已报道的含能分子中获取稠环类分子骨架数据集，然后从中筛选出具有高密度潜质的分子骨架与氨基和硝基进行片段组装，最终得到大量的待筛选分子结构组成虚拟筛选空间。在性能预测阶段，首先采用图神经网络训练的QSPR 模型预测晶体密度，然后针对晶体密度高的分子采用高通量计算平台EM-Studio 预测生成焓、BDE、爆速、爆压和撞击感度（H50），并由综合性能筛选获得含能分子候选物。

1.2 研究方法

1.2.1 虚拟筛选空间构建

（1）分子骨架获取

通过收集文献报道的元素范围为CHON 的硝基化合物为基本数据集，将所有分子的基团逐一摘除后获得有机分子骨架（图2）。

图2 获取分子骨架的原理Fig.2 Principle for getting molecular skeletons

（2）高密度分子骨架筛选

分子骨架的密度采用分子等效密度来评估，其中，等效密度的计算通过式（1）得到：

式中，M为相对分子质量，g·mol-1；V为分子体积，A3，其中1 Mol分子中含有Na（阿伏伽德罗常数，6.022×1023）个分子，并且需要将体积换算成立方厘米（cm3，1 cm3=1×1024A3），经过单位换算，得到

基于式（2），筛选出密度大于1.8 g·cm-3的有机骨架用于虚拟筛选空间的构建。

（3）分子片段组装

为了让分子结构可能具有含能分子的性质，需要将分子骨架与致爆基团和稳定基团进行片段组装。研究中以硝基为致爆基团，以氨基为稳定基团，组装方式为分子骨架的不同位点被氨基或硝基取代。分子片段组装采用了深度优先搜索算法，该方法的本质是用分子的所有可取代位点构建一颗搜索树，通过对搜索树树枝的遍历，获取出需要取代的位点进行取代。

（4）稠环结构的识别

研究表明和传统含能分子相比，稠环含能分子因具有较高的生成热和较大的环张力，使得该类分子具有较为优异的爆轰性能［35-36］。因此在选择分子骨架的时候，应尽可能的寻找稠环分子骨架，保证设计出来的分子为稠环含能分子。本研究采用的是Tarjan 算法［37］来识别稠环骨架，如图3 所示。

如图3 所示，Tarjan 算法将分子处理为一个图结构，以原子作为图的顶点，以化学键作为图的边。对每个点进行搜索序和回溯值的计算，最终确定出图中的环的数量（图中红线的数量），将并环数量大于2 的结构识别为稠环。将左边的分子结构进行编号，通过深度优先搜索算法对其进行遍历，遍历后得到搜索树，对每个节点进行回溯值进行计算，就可以确定环的数量和大小。

图3 Tarjan 算法求环的过程Fig.3 Tarjan algorithm to find a ring

1.2.2 性能预测

（1）密度的预测

密度是影响含能分子爆轰性能的一个关键属性，准确的密度预测模型对于含能分子最终的筛选成功率将起到至关重要的影响。近年来，很多机器学习和深度学习的方法在预测分子性质中得到运用，本研究用于含能分子密度预测的深度学习模型采用了基于MEGNet［38］的图神经网络模型训练密度预测模型。该方法于2018 年提出，可以通过分子结构直接预测出分子的性质，性能较好，该模型包含特征提取，特征向量化和模型回归3 个步骤（图4）。首先将分子结构读入网络中，用RDKit 获取分子的特征矩阵，包括3D 距离矩阵等信息，将这些信息通过图卷积网络得到特征向量，最后通过多层感知机来计算分子密度。为了验证模型预测的准确性，需要计算模型预测的精度，本研究将数据集分成80%的训练集，10%的验证集和10%的独立验证集，模型的误差分析采用平均绝对误差（MAE）、均方根误差（RMSE）和决定系数（R2），计算过程如公式（3）～（6）所示。

图4 基于MEGNet 的密度预测模型框架Fig.4 Framework of density prediction model based on MEGNet

（4）硝基键BDE 计算

根据“引发键”理论［40-42］，分子骨架与硝基直接相连的化学键是含能分子分解反应最容易先发生断裂的键，因此该硝基键的BDE 常用来评估分子稳定性。研究计算BDE 采用的公式为：

式中，BDE为分子的键解离能，kJ·mol-1，E（R—NO2）为分子本身的能量，kJ·mol-1，E(·NO2)为硝基自由基的能量，kJ·mol-1，E(·R)为分子减去硝基后的自由基的能量，kJ·mol-1。

（5）撞击感度预测

含能化合物的撞击感度也是常见的稳定性的判据，采用Keshavarz［45］的方法计算撞击感度：

式中，H50为特性落高［43］，a为分子中C 原子的数量，b为分子中H 原子的数量，c为分子中N 原子的数量，d为分子中O 原子的数量，M为分子的相对分子质量。

1.3 软件与工具

研究中的GNN 模型基于Chainer 架构构建，使用与MEGNet［38］相同的图块，在训练过程中采用Adam optimizer 和ReLu［44］激 活 函 数。当 前DL 研 究 的 运 行环 境 是Python 虚 拟 环 境，嵌 入 了RDkit［45］和Chainer_Chemistry［46］。密度预测模型的数据集从CCDC 数据库［47］中采集得到，总计超过2000 个硝基化合物。所有量子化学计算由Gaussian09 程序执行，量子化学计算的前后的处理、生成焓、BDE 爆速和爆压的计算由高通量计算平台EM-Studio 自动完成。

表1 K-J 方程参数的计算Table1 Calculation of the parameters of K-J equation

2 结果与讨论

2.1 虚拟筛选搜索空间构建

研究设想通过预筛选高密度骨架为后续获得高密度的待筛选分子对象，因此预先分析了已报道分子的骨架等效密度与对应含能分子的晶体密度之间的关系（图5）。图5 表明，含能分子的晶体密度随其骨架密度的增加而增加的大致趋势，两者的线性相关度R达到了0.4，表明含能分子的晶体密度与其骨架密度具有中度的正相关性。由此可知，基于高密度骨架来构建具有潜在高密度的待筛选分子可行。

图5 分子骨架的等效密度与晶体密度之间的关系，横坐标为分子骨架的等效密度，纵坐标为分子的晶体密度Fig.5 Relationship between the equivalent density of molecular skeleton and the crystal density.The abscissa is the equivalent density of molecular skeleton（ρske）and the ordinate is the crystal density of molecule（ρcrst）

在虚拟筛选空间的构建中，为了避免设计出来的分子都集中在某一种或某几种结构中，需要保证分子结构的多样性，以实现设计出来的分子能够覆盖大多数结构类型。研究以分子骨架选择的多样性作为实现分子的多样性的前提，基于分子密度、原子总数、氮元素的含量和分子中环的个数这4 个指标筛选出等效密度在1.8～2.4 之间的分子骨架总计200 个（图6）。

由图6 可以看出，骨架分子包含的原子总数范围为10～22 个之间，并且原子数量集中在17～21 个，整体原子数量符合正态分布的规律，说明筛选出的分子骨架符合多样性的特点，在骨架大小上具有明显的差异。分子中氮含量的质量百分比范围为8%～80%，说明筛选出的分子骨架在氮含量上具有明显差异，氮含量的差异可能会直接影响设计出来的分子爆轰性能。环数集中在2～4 个之间，多个环是稠环的基础，只有保证分子中含有多个环，才能保证最终设计的分子有可能为稠环分子。通过图6 我们看到，现有分子的骨架符合分子多样性特点。研究基于上述200 个骨架结构，以硝基为致爆基团，以氨基为稳定基团，通过分子组装获得了近7000 个分子构成虚拟筛选空间。

图6 虚拟筛选空间分析Fig.6 Analysis of the virtual screening space

2.2 密度预测模型

本研究为了设计数据集，从CCDC 中搜索到含有硝基的CHON 类分子，然后选出具有常温常压下有密度数据的分子总计2000，作为密度预测模型训练的数据集。将整个数据集分别分为80%的训练集、10% 的验证集和10% 的独立验证集，对于GNN 模型，只有训练集参与训练，验证集用于判断收敛点和防止过拟合。为了保证独立验证集的独立性，对独立验证集进行了Tanimoto 相似性［48］分析，Tamnimoto相似性是将分子转换成二进制数，判断2 个分子中的二进制数相同的个数，这种方法可以有效的区分相似分子，相同的独立测试集中的样本与其在训练/验证集中的最近邻之间的Tanimoto 相似性的统计结果如图7 所示。

图7 验证集与其训练集中的最近邻间的相似性，Ej，ave和Ej，std分别代表Ej的平均值和标准偏差Fig.7 Similarity between the samples in the test set and their nearest neighbors in the training/verification set. Ej，ave and Ej，std represent the mean and standard deviation of Ej，respectively

如图7 所示，测试集中每个测试样本与其训练/验证集中的最近邻的Tanimoto 相似性Ej在统计上呈正态分布，标准偏差Ej，std等于0.069，表明样本分布符合随机抽样的特征。平均Tanimoto 相似性Ej，ave=0.71，表明测试集与训练/验证集具有相关性的同时还具有相对的独立性。图7 的结果表明，密度预测模型所采用的集试是独立于训练/验证集的，可满足模型测试对独立测试集的要求，为确保模型的泛化性提供了前提。

基于上述结果，我们将测试集中的样本用于由MEGNET 训练得到的密度预测模型，测试其精度与泛化性（图8）。

图8 （a）密度的预测结果宇称图和（b）误差分布图Fig.8 （a）Parity diagram of density prediction results and（b）error distribution diagram

如 图8a 所 示，模 型 的MAE、RMSE 和R2分 别 为0.044、0.059 和0.924。如图8b 所示，模型的预测误差（Δρ）主要分布在-5%和5%之间，误差分布集中在0%～3%附近。这一结果表明，本研究得到的基于图神经网络的密度预测模型具有可靠的精度，能够用于后续以晶体密度为指标的含能分子筛选。

2.3 高密度含能分子的筛选

为了评价含能分子的筛选效果，这里将采用密度的指标来进行评价，结果如图9 所示。

图9 虚拟筛选空间样本的密度分布图Fig.9 Density distribution of samples in the virtual screening space

如图9 所示可以看到，筛选空间中的分子的预测晶体密度全部高于1.65 g·cm-3，且预测晶体密处于1.8～1.9 g·cm-3范围的分子占到了大多数，晶体密度接近2.0 g·cm-3的分子数量达到了28。这一结果表明，采用预筛选高密度骨架的方式构建的含能分子筛选空间，其分子样本具有显著的高晶体密度特征。这一结果验证了图5 中的分析结果，即分子骨架等效密度与晶体密度具有正相关性，预筛选高密度分子骨架可以有效提升高密度含能分子的筛选成功率。

为了评估本研究设计的含能分子的综合性能，选取密度排名前100 的分子，由EM-Studio 计算其生成焓、爆速、爆压和BDE。为了直观体现被筛选分子的综合性能，采用常见的含能分子TNT、RDX 和TATB 作为指标分子，对各种性质进行了对比分析（图10）。

图10 候选含能分子的性能，（a）爆压，（b）爆速，（c）BDE 与（d）晶体密度（红色线为TATB 的对应性质，绿色线为TNT 的对应性质，蓝色线为RDX 的对应性质）。Fig.10 Properties of the candidate energetic molecules，（a）detonation pressure，（b）detonation velocity，（c）BDE，and（d）crystal density. The red line is the corresponding property of TATB，the green line is the corresponding property of TNT，and the blue line is the corresponding property of RDX

从图10 中我们可以看到，有超过半数的分子爆压超过TATB（图a），并且爆速整体超过TATB，有部分分子爆速甚至超过RDX（图b），在BDE 上，大部分分子超过RDX，过半分子超过TNT（图c），全部分子密度超过TNT，并且大部分超过RDX（图d），这说明通过密度筛选出来的分子的爆轰能力普遍大于TATB 和TNT，并且个别分子甚至超过了RDX，可能出现的原因是因为实验筛选出的分子的氮含量都比较高，且环的数量都大于等于2，这类分子被认为［49］环内的N—N 键和C—N 键在爆炸时会释放出大量的能量，可有效提升爆轰性能。从图10c 中可以看到有接近一半的分子稳定性优于RDX，其原因可能是因为分子骨架中大量含有呋咱、四嗪、四唑和三唑等富氮杂环，研究表明［50］，这些分子具有良好的热稳定性，已经成为部分研究设计含能分子的有效结构单元。

在筛选分子的时候将爆速以RDX 为基准，BDE 以TNT 为基准进行对比，以期得到爆轰性能强于RDX，且稳定性强于TNT 的含能分子（图11）。

如图11 所示，在BDE 与爆速同时较高的分子数量较少，爆速超过RDX，BDE 超过TNT 的分子仅仅只有6 个，进一步验证了含能分子的能量与安全性之间可能存在一些矛盾，高能量的分子往往稳定性较差，在设计分子的时候应高要综合考虑爆轰性能与稳定性的关系。研究表明［3］，在分子水平上，能量和稳定性的矛盾尤为突出，但在晶体和混合物层次，能量和安全性的矛盾可以得到很大程度上的缓解。因此随着后期的晶体工程和复合工艺的发展，可能有希望得到越来越多的高能量高稳定性的含能分子。这6个分子的环中的化学键均为大π 键，并且多个环之间仅有一条化学键相连，或许这样的结构能在保证分子化学稳定性的基础之上，较大的提升分子在爆轰时释放出来的能量，后续在设计分子和选择骨架的时候，应多选择此类骨架和分子。具体分子性质如表2 所示。

图11 Ⅰ区为候选分子爆速与BDE 的权衡，Ⅱ区为在Ⅰ区中右上角爆速大于RDX 且BDE 大于TNT 的分子结构Fig.11 Panel Ⅰis the tradeoff between candidate molecular detonation velocity and BDE. Panel Ⅱis the molecular structure in which the detonation velocity in the upper right corner of Panel Ⅰis greater than RDX and BDE is greater than TNT

表2 中可以看到生成的分子在爆轰性质上基本超过TATB，在部分性质上甚至超过RDX，化学稳定性均强于RDX，因此可以尝试将本次实验设计的化合物进行进一步的研究。值得关注的是，本研究设计的分子在稳定性上仍无法超过TATB。TATB 在结构上具有独有的特点，其3 个硝基和3 个氨基是相间排列在苯环上的，使整个分子体系的电荷分布高度对称，相间分布的氨基和硝基通过强的分子内氢键和静电作用在分子外围形成非常稳定的环状结构［53］。从当前的研究来看，设计出类似于TATB 这样完美的分子结构非常困难，这可能是探索新型高稳定性含能分子需要面对的挑战。

表2 含能分子候选物的性能Table 2 Properties of the candidates for energetic molecules

3 结论

（1）含能分子的晶体密度与骨架密度成正相关，基于高密度骨架可构建具有潜在高密度的待筛选分子。

（2）通过分子骨架预筛选获得200 个高密度的稠环分子骨架，用硝基和氨基作为取代基，可以组合出7000 多个潜在的含能分子，且均具备高密度的潜质。

（3）基于CCDC 数据库中的硝基化合物为数据集，通过图神经网络训练出晶体密度预测模型，其决定系数R2达到0.924，且泛化性良好。

（4）基于深度学习、量子化学计算和爆轰产物状态方程等方法预测晶体密度、生成焓、爆轰性能和化学稳定性，从而由性能排序筛选出能量水平优于RDX，稳定性优于TNT 的新型含能分子6 个。

（5）通过预筛选分子骨架可以有效提升虚拟筛选空间的总体性能，在此基础上可借助高通量计算与深度学习实现含能分子的高效设计。