王潤文，楊春明，劉 建

（1. 西南科技大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，四川 綿陽 621000；2. 中國工程物理研究院化工材料研究所，四川 綿陽 621999）

0 引言

含能分子是一類具有高能量屬性的特殊分子，廣泛用于國防、航空航天和民用領(lǐng)域［1-4］。含能分子在研制的各環(huán)節(jié)中必須兼顧多項(xiàng)性能指標(biāo)，包括能量［5-7］、吸濕性［8-10］、熱穩(wěn)定性［11-12］、化學(xué)鍵能［13-14］、感度［15-16］等。因此，新含能分子從最初的設(shè)計(jì)到最終的應(yīng)用需要經(jīng)歷復(fù)雜的流程，發(fā)現(xiàn)真正有應(yīng)用前景的新含能分子極具挑戰(zhàn)性。20 世紀(jì)80 年代開始，計(jì)算機(jī)輔助材料設(shè)計(jì)（Compute-aided material design，CAMD）［17］開始興起，并逐步被引入到含能材料的研究中，使含能材料的性能可以通過計(jì)算機(jī)手段來獲得［18-21］。CAMD 可以通過預(yù)測(cè)含能分子的性質(zhì)提前評(píng)估其是否具備合成價(jià)值，從而篩選出有價(jià)值的分子，減少了實(shí)驗(yàn)中的“試錯(cuò)”行為。在CAMD 的協(xié)助下，含能分子的研制模式已由傳統(tǒng)的“試錯(cuò)法”升級(jí)為計(jì)算與實(shí)驗(yàn)協(xié)同的模式，即先由計(jì)算篩選出滿足性質(zhì)的目標(biāo)含能分子，然后探索合成路線。得益于此，化學(xué)家可以避免在無價(jià)值的分子上耗費(fèi)時(shí)間和資源，使含能分子的研制效率在一定程度上得到提升。

在實(shí)際應(yīng)用中，含能分子設(shè)計(jì)的成功率取決于2個(gè)主要因素，即分子結(jié)構(gòu)的構(gòu)建方式和性能預(yù)測(cè)方法。在分子結(jié)構(gòu)的構(gòu)建方面，由于不同的元素組成和原子間連接存在多種可能性，CAMD 的研究對(duì)象往往是由一定數(shù)量的分子結(jié)構(gòu)組成集合，也即虛擬篩選空間［22］。在含能分子的性能預(yù)測(cè)方面，一般采用量子化學(xué)計(jì)算或定量構(gòu)效關(guān)系（Quantitative structure-property relationship，QSPR）模型［19，23-24］來實(shí)現(xiàn)。含能分子設(shè)計(jì)的本質(zhì)就是基于性能預(yù)測(cè)的結(jié)果對(duì)虛擬篩選空間中的分子結(jié)構(gòu)進(jìn)行排序，排序靠前的分子結(jié)構(gòu)將作為候選含能分子推薦給合成化學(xué)家。因此，為保證含能分子設(shè)計(jì)的“成功率”，高效的虛擬篩選空間構(gòu)建方式和準(zhǔn)確的性能預(yù)測(cè)方法是基本前提。近年來，基于機(jī)器學(xué)習(xí)的QSPR 和高通量篩選技術(shù)開始用于CAMD，并在含能分子設(shè)計(jì)中獲得了廣泛的應(yīng)用，使含能分子的設(shè)計(jì)水平大幅提升。

Pinheiro 等［25］通過機(jī)器學(xué)習(xí)預(yù)測(cè)QM9 數(shù)據(jù) 集中分子的9 種性質(zhì)，最終模型精度最高HUMO-LUMO預(yù)測(cè)精度（R2）達(dá)到了0.97。Kang 等［26］通過機(jī)器學(xué)習(xí)預(yù)測(cè)含能分子的生成焓等性質(zhì)，其訓(xùn)練的模型的R2達(dá)到0.93。Elton 等［27］以鍵總數(shù)為描述符，采用Kernel Ridge 回歸，得到晶體密度性質(zhì)的預(yù)測(cè)模型，R2為0.74。Casey 等［28］采用分子三維電子結(jié)構(gòu)為描述符，基于Deep Neural Network（DNN）得到含能材料7 種主要性質(zhì)的預(yù)測(cè)模型，7 種模型的R2值均達(dá)到0.95 以上。徐雅斌等［29］以庫侖矩陣為輸入，基于Convolutional Neural Network（CNN）和Bi-directional Short-term Memory Network（Bi-LSTM）的融合模型，實(shí)現(xiàn)了對(duì)含能材料生成焓的預(yù)測(cè)，其訓(xùn)練出的模型的R2為0.97。Yang 等［30］基于圖神經(jīng)網(wǎng)絡(luò)訓(xùn)練了一個(gè)密度預(yù)測(cè)模型，可以通過分子結(jié)構(gòu)直接預(yù)測(cè)出晶體密度，模型精度R2值達(dá)到了0.949。Jiang 等［31］通過神經(jīng)網(wǎng)絡(luò)獲得了分子的晶體的共晶可行性判據(jù)模型，在此基礎(chǔ)上成功制備出一種新型含能共晶。Huang 等［32］采用晶體描述符和梯度提升機(jī)（XGBoost）等5 種算法獲得了面向爆熱、爆速、爆壓、熱分解溫度和晶格能的機(jī)器學(xué)習(xí)模型，建立了綜合評(píng)估爆轟性能和穩(wěn)定性的高通量解決方案。Liu 等［33］運(yùn)用高通量篩選技術(shù)設(shè)計(jì)并篩選了大量含 有CHONF 的 含 能 分 子。Song 等［34］基 于 機(jī) 器 學(xué) 習(xí)建立了分子性質(zhì)預(yù)測(cè)和晶體堆積模式評(píng)估模型，從25 112 個(gè)分子構(gòu)成的虛擬篩選空間中篩選出綜合性優(yōu)異的新型含能分子，并實(shí)現(xiàn)了實(shí)驗(yàn)室合成與測(cè)試。

上述研究工作將機(jī)器學(xué)習(xí)和高通量篩選技術(shù)用于含能分子的設(shè)計(jì)，為含能分子的高效篩選提供了解決方案，但在虛擬篩選空間構(gòu)建和性能預(yù)測(cè)模型方面仍有改進(jìn)的空間。對(duì)于虛擬篩選空間構(gòu)建，當(dāng)前的CAMD 需要研究者依據(jù)知識(shí)和經(jīng)驗(yàn)選定基本結(jié)構(gòu)單元和官能團(tuán)種類后再逐一枚舉分子結(jié)構(gòu)。該方法以人為判斷為前提，易造成虛擬篩選空間構(gòu)建效率低和分子結(jié)構(gòu)類型覆蓋面窄等問題。對(duì)于性能預(yù)測(cè)方法，有多種方法可以實(shí)現(xiàn)對(duì)爆轟性能準(zhǔn)確的預(yù)測(cè)，且與爆轟性能密切相關(guān)的密度和生成焓預(yù)測(cè)結(jié)果一般較可靠。存在的問題是部分模型基于的訓(xùn)練樣本少、適用范圍窄，尤其是化學(xué)穩(wěn)定性相關(guān)的性質(zhì)（比如感度）一直以來都難以被準(zhǔn)確預(yù)測(cè)。

本研究擬將含能分子研究中廣受關(guān)注的稠環(huán)類結(jié)構(gòu)為研究對(duì)象，開展高通量計(jì)算與深度學(xué)習(xí)相結(jié)合的含能分子設(shè)計(jì)，并針對(duì)上述問題提出相應(yīng)的解決方案。在虛擬篩選空間構(gòu)建中，基于分子骨架等效密度的預(yù)篩選獲得了高密度分子骨架，然后再通過分子片段組裝構(gòu)建了具有潛在高密度特征的虛擬篩選空間。在性能預(yù)測(cè)中，根據(jù)不同預(yù)測(cè)方法的適用性，（i）采用量子化學(xué)計(jì)算預(yù)測(cè)生成焓、X-NO2鍵解離能（Bond Dissociation Energy，BDE），（ii）采用深度學(xué)習(xí)模型預(yù)測(cè)晶體密度；（iii）采用爆轟性能經(jīng)驗(yàn)公式預(yù)測(cè)爆轟性能。針對(duì)含能分子化學(xué)穩(wěn)定性難以準(zhǔn)確預(yù)測(cè)的問題，采用與化學(xué)穩(wěn)定性密切相關(guān)的X-NO2BDE［56］間接評(píng)估化學(xué)穩(wěn)定性。在高密度骨架預(yù)篩選的基礎(chǔ)上構(gòu)建了具有潛在高密度的虛擬篩選空間，并將高通量計(jì)算和機(jī)器學(xué)習(xí)相結(jié)合預(yù)測(cè)了密度、生成焓、爆速、爆壓和X-NO2BDE，從中篩選出性能優(yōu)異的候選含能分子共計(jì)6 種。

1 原理與方法

1.1 基本原理

基于CAMD 的含能分子篩選，首先需要構(gòu)建由大量待篩選分子構(gòu)成的虛擬篩選空間，然后采用量子化學(xué)計(jì)算或是QSPR 模型進(jìn)行篩選，最終獲得符合性能預(yù)期的含能分子。高效的分子篩選需要解決2 個(gè)問題：一是如何構(gòu)建具有高性能潛質(zhì)篩選對(duì)象，以提高分子篩選的成功率；二是如何獲得精度高、泛化性好的模型以實(shí)現(xiàn)對(duì)性能的準(zhǔn)確預(yù)測(cè)。本研究以近年來較熱門的稠環(huán)類分子作為研究對(duì)象，基于高通量計(jì)算與深度學(xué)習(xí)相結(jié)合的方法，力圖發(fā)現(xiàn)高密度、高爆轟性能和高化學(xué)穩(wěn)定性的含能分子候選物。具體的技術(shù)原理如圖1 所示。

圖1 高通量計(jì)算與深度學(xué)習(xí)相結(jié)合的稠環(huán)含能分子設(shè)計(jì)技術(shù)原理Fig.1 Principle for designing fused-ring energetic compounds using high-throughput computing combined with deep learning

如圖1 所示，研究首先預(yù)先篩選高密度分子骨架，然后經(jīng)分子片段組裝獲得由潛在的高密度含能分子構(gòu)成的虛擬篩選空間。在此基礎(chǔ)上基于深度學(xué)習(xí)、量子化學(xué)計(jì)算和爆轟產(chǎn)物狀態(tài)方程等方法預(yù)測(cè)晶體密度、生成焓、爆轟性能、BDE 和撞擊感度的數(shù)值，從而由性能排序完成含能分子的篩選。虛擬篩選空間的構(gòu)建由分子高密度骨架篩選和分子片段對(duì)接兩步構(gòu)成。首先，從收集到的已報(bào)道的含能分子中獲取稠環(huán)類分子骨架數(shù)據(jù)集，然后從中篩選出具有高密度潛質(zhì)的分子骨架與氨基和硝基進(jìn)行片段組裝，最終得到大量的待篩選分子結(jié)構(gòu)組成虛擬篩選空間。在性能預(yù)測(cè)階段，首先采用圖神經(jīng)網(wǎng)絡(luò)訓(xùn)練的QSPR 模型預(yù)測(cè)晶體密度，然后針對(duì)晶體密度高的分子采用高通量計(jì)算平臺(tái)EM-Studio 預(yù)測(cè)生成焓、BDE、爆速、爆壓和撞擊感度（H50），并由綜合性能篩選獲得含能分子候選物。

1.2 研究方法

1.2.1 虛擬篩選空間構(gòu)建

（1）分子骨架獲取

通過收集文獻(xiàn)報(bào)道的元素范圍為CHON 的硝基化合物為基本數(shù)據(jù)集，將所有分子的基團(tuán)逐一摘除后獲得有機(jī)分子骨架（圖2）。

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

（2）高密度分子骨架篩選

分子骨架的密度采用分子等效密度來評(píng)估，其中，等效密度的計(jì)算通過式（1）得到：

式中，M為相對(duì)分子質(zhì)量，g·mol-1；V為分子體積，A3，其中1 Mol分子中含有Na（阿伏伽德羅常數(shù)，6.022×1023）個(gè)分子，并且需要將體積換算成立方厘米（cm3，1 cm3=1×1024A3），經(jīng)過單位換算，得到

基于式（2），篩選出密度大于1.8 g·cm-3的有機(jī)骨架用于虛擬篩選空間的構(gòu)建。

（3）分子片段組裝

為了讓分子結(jié)構(gòu)可能具有含能分子的性質(zhì)，需要將分子骨架與致爆基團(tuán)和穩(wěn)定基團(tuán)進(jìn)行片段組裝。研究中以硝基為致爆基團(tuán)，以氨基為穩(wěn)定基團(tuán)，組裝方式為分子骨架的不同位點(diǎn)被氨基或硝基取代。分子片段組裝采用了深度優(yōu)先搜索算法，該方法的本質(zhì)是用分子的所有可取代位點(diǎn)構(gòu)建一顆搜索樹，通過對(duì)搜索樹樹枝的遍歷，獲取出需要取代的位點(diǎn)進(jìn)行取代。

（4）稠環(huán)結(jié)構(gòu)的識(shí)別

研究表明和傳統(tǒng)含能分子相比，稠環(huán)含能分子因具有較高的生成熱和較大的環(huán)張力，使得該類分子具有較為優(yōu)異的爆轟性能［35-36］。因此在選擇分子骨架的時(shí)候，應(yīng)盡可能的尋找稠環(huán)分子骨架，保證設(shè)計(jì)出來的分子為稠環(huán)含能分子。本研究采用的是Tarjan 算法［37］來識(shí)別稠環(huán)骨架，如圖3 所示。

如圖3 所示，Tarjan 算法將分子處理為一個(gè)圖結(jié)構(gòu)，以原子作為圖的頂點(diǎn)，以化學(xué)鍵作為圖的邊。對(duì)每個(gè)點(diǎn)進(jìn)行搜索序和回溯值的計(jì)算，最終確定出圖中的環(huán)的數(shù)量（圖中紅線的數(shù)量），將并環(huán)數(shù)量大于2 的結(jié)構(gòu)識(shí)別為稠環(huán)。將左邊的分子結(jié)構(gòu)進(jìn)行編號(hào)，通過深度優(yōu)先搜索算法對(duì)其進(jìn)行遍歷，遍歷后得到搜索樹，對(duì)每個(gè)節(jié)點(diǎn)進(jìn)行回溯值進(jìn)行計(jì)算，就可以確定環(huán)的數(shù)量和大小。

圖3 Tarjan 算法求環(huán)的過程Fig.3 Tarjan algorithm to find a ring

1.2.2 性能預(yù)測(cè)

（1）密度的預(yù)測(cè)

密度是影響含能分子爆轟性能的一個(gè)關(guān)鍵屬性，準(zhǔn)確的密度預(yù)測(cè)模型對(duì)于含能分子最終的篩選成功率將起到至關(guān)重要的影響。近年來，很多機(jī)器學(xué)習(xí)和深度學(xué)習(xí)的方法在預(yù)測(cè)分子性質(zhì)中得到運(yùn)用，本研究用于含能分子密度預(yù)測(cè)的深度學(xué)習(xí)模型采用了基于MEGNet［38］的圖神經(jīng)網(wǎng)絡(luò)模型訓(xùn)練密度預(yù)測(cè)模型。該方法于2018 年提出，可以通過分子結(jié)構(gòu)直接預(yù)測(cè)出分子的性質(zhì)，性能較好，該模型包含特征提取，特征向量化和模型回歸3 個(gè)步驟（圖4）。首先將分子結(jié)構(gòu)讀入網(wǎng)絡(luò)中，用RDKit 獲取分子的特征矩陣，包括3D 距離矩陣等信息，將這些信息通過圖卷積網(wǎng)絡(luò)得到特征向量，最后通過多層感知機(jī)來計(jì)算分子密度。為了驗(yàn)證模型預(yù)測(cè)的準(zhǔn)確性，需要計(jì)算模型預(yù)測(cè)的精度，本研究將數(shù)據(jù)集分成80%的訓(xùn)練集，10%的驗(yàn)證集和10%的獨(dú)立驗(yàn)證集，模型的誤差分析采用平均絕對(duì)誤差（MAE）、均方根誤差（RMSE）和決定系數(shù)（R2），計(jì)算過程如公式（3）～（6）所示。

圖4 基于MEGNet 的密度預(yù)測(cè)模型框架Fig.4 Framework of density prediction model based on MEGNet

（4）硝基鍵BDE 計(jì)算

根據(jù)“引發(fā)鍵”理論［40-42］，分子骨架與硝基直接相連的化學(xué)鍵是含能分子分解反應(yīng)最容易先發(fā)生斷裂的鍵，因此該硝基鍵的BDE 常用來評(píng)估分子穩(wěn)定性。研究計(jì)算BDE 采用的公式為：

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

（5）撞擊感度預(yù)測(cè)

含能化合物的撞擊感度也是常見的穩(wěn)定性的判據(jù)，采用Keshavarz［45］的方法計(jì)算撞擊感度：

式中，H50為特性落高［43］，a為分子中C 原子的數(shù)量，b為分子中H 原子的數(shù)量，c為分子中N 原子的數(shù)量，d為分子中O 原子的數(shù)量，M為分子的相對(duì)分子質(zhì)量。

1.3 軟件與工具

研究中的GNN 模型基于Chainer 架構(gòu)構(gòu)建，使用與MEGNet［38］相同的圖塊，在訓(xùn)練過程中采用Adam optimizer 和ReLu［44］激 活 函 數(shù)。當(dāng) 前DL 研 究 的 運(yùn) 行環(huán) 境 是Python 虛 擬 環(huán) 境，嵌 入 了RDkit［45］和Chainer_Chemistry［46］。密度預(yù)測(cè)模型的數(shù)據(jù)集從CCDC 數(shù)據(jù)庫［47］中采集得到，總計(jì)超過2000 個(gè)硝基化合物。所有量子化學(xué)計(jì)算由Gaussian09 程序執(zhí)行，量子化學(xué)計(jì)算的前后的處理、生成焓、BDE 爆速和爆壓的計(jì)算由高通量計(jì)算平臺(tái)EM-Studio 自動(dòng)完成。

表1 K-J 方程參數(shù)的計(jì)算Table1 Calculation of the parameters of K-J equation

2 結(jié)果與討論

2.1 虛擬篩選搜索空間構(gòu)建

研究設(shè)想通過預(yù)篩選高密度骨架為后續(xù)獲得高密度的待篩選分子對(duì)象，因此預(yù)先分析了已報(bào)道分子的骨架等效密度與對(duì)應(yīng)含能分子的晶體密度之間的關(guān)系（圖5）。圖5 表明，含能分子的晶體密度隨其骨架密度的增加而增加的大致趨勢(shì)，兩者的線性相關(guān)度R達(dá)到了0.4，表明含能分子的晶體密度與其骨架密度具有中度的正相關(guān)性。由此可知，基于高密度骨架來構(gòu)建具有潛在高密度的待篩選分子可行。

圖5 分子骨架的等效密度與晶體密度之間的關(guān)系，橫坐標(biāo)為分子骨架的等效密度，縱坐標(biāo)為分子的晶體密度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）

在虛擬篩選空間的構(gòu)建中，為了避免設(shè)計(jì)出來的分子都集中在某一種或某幾種結(jié)構(gòu)中，需要保證分子結(jié)構(gòu)的多樣性，以實(shí)現(xiàn)設(shè)計(jì)出來的分子能夠覆蓋大多數(shù)結(jié)構(gòu)類型。研究以分子骨架選擇的多樣性作為實(shí)現(xiàn)分子的多樣性的前提，基于分子密度、原子總數(shù)、氮元素的含量和分子中環(huán)的個(gè)數(shù)這4 個(gè)指標(biāo)篩選出等效密度在1.8～2.4 之間的分子骨架總計(jì)200 個(gè)（圖6）。

由圖6 可以看出，骨架分子包含的原子總數(shù)范圍為10～22 個(gè)之間，并且原子數(shù)量集中在17～21 個(gè)，整體原子數(shù)量符合正態(tài)分布的規(guī)律，說明篩選出的分子骨架符合多樣性的特點(diǎn)，在骨架大小上具有明顯的差異。分子中氮含量的質(zhì)量百分比范圍為8%～80%，說明篩選出的分子骨架在氮含量上具有明顯差異，氮含量的差異可能會(huì)直接影響設(shè)計(jì)出來的分子爆轟性能。環(huán)數(shù)集中在2～4 個(gè)之間，多個(gè)環(huán)是稠環(huán)的基礎(chǔ)，只有保證分子中含有多個(gè)環(huán)，才能保證最終設(shè)計(jì)的分子有可能為稠環(huán)分子。通過圖6 我們看到，現(xiàn)有分子的骨架符合分子多樣性特點(diǎn)。研究基于上述200 個(gè)骨架結(jié)構(gòu)，以硝基為致爆基團(tuán)，以氨基為穩(wěn)定基團(tuán)，通過分子組裝獲得了近7000 個(gè)分子構(gòu)成虛擬篩選空間。

圖6 虛擬篩選空間分析Fig.6 Analysis of the virtual screening space

2.2 密度預(yù)測(cè)模型

本研究為了設(shè)計(jì)數(shù)據(jù)集，從CCDC 中搜索到含有硝基的CHON 類分子，然后選出具有常溫常壓下有密度數(shù)據(jù)的分子總計(jì)2000，作為密度預(yù)測(cè)模型訓(xùn)練的數(shù)據(jù)集。將整個(gè)數(shù)據(jù)集分別分為80%的訓(xùn)練集、10% 的驗(yàn)證集和10% 的獨(dú)立驗(yàn)證集，對(duì)于GNN 模型，只有訓(xùn)練集參與訓(xùn)練，驗(yàn)證集用于判斷收斂點(diǎn)和防止過擬合。為了保證獨(dú)立驗(yàn)證集的獨(dú)立性，對(duì)獨(dú)立驗(yàn)證集進(jìn)行了Tanimoto 相似性［48］分析，Tamnimoto相似性是將分子轉(zhuǎn)換成二進(jìn)制數(shù)，判斷2 個(gè)分子中的二進(jìn)制數(shù)相同的個(gè)數(shù)，這種方法可以有效的區(qū)分相似分子，相同的獨(dú)立測(cè)試集中的樣本與其在訓(xùn)練/驗(yàn)證集中的最近鄰之間的Tanimoto 相似性的統(tǒng)計(jì)結(jié)果如圖7 所示。

圖7 驗(yàn)證集與其訓(xùn)練集中的最近鄰間的相似性，Ej，ave和Ej，std分別代表Ej的平均值和標(biāo)準(zhǔn)偏差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 所示，測(cè)試集中每個(gè)測(cè)試樣本與其訓(xùn)練/驗(yàn)證集中的最近鄰的Tanimoto 相似性Ej在統(tǒng)計(jì)上呈正態(tài)分布，標(biāo)準(zhǔn)偏差Ej，std等于0.069，表明樣本分布符合隨機(jī)抽樣的特征。平均Tanimoto 相似性Ej，ave=0.71，表明測(cè)試集與訓(xùn)練/驗(yàn)證集具有相關(guān)性的同時(shí)還具有相對(duì)的獨(dú)立性。圖7 的結(jié)果表明，密度預(yù)測(cè)模型所采用的集試是獨(dú)立于訓(xùn)練/驗(yàn)證集的，可滿足模型測(cè)試對(duì)獨(dú)立測(cè)試集的要求，為確保模型的泛化性提供了前提。

基于上述結(jié)果，我們將測(cè)試集中的樣本用于由MEGNET 訓(xùn)練得到的密度預(yù)測(cè)模型，測(cè)試其精度與泛化性（圖8）。

圖8 （a）密度的預(yù)測(cè)結(jié)果宇稱圖和（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 所示，模型的預(yù)測(cè)誤差（Δρ）主要分布在-5%和5%之間，誤差分布集中在0%～3%附近。這一結(jié)果表明，本研究得到的基于圖神經(jīng)網(wǎng)絡(luò)的密度預(yù)測(cè)模型具有可靠的精度，能夠用于后續(xù)以晶體密度為指標(biāo)的含能分子篩選。

2.3 高密度含能分子的篩選

為了評(píng)價(jià)含能分子的篩選效果，這里將采用密度的指標(biāo)來進(jìn)行評(píng)價(jià)，結(jié)果如圖9 所示。

圖9 虛擬篩選空間樣本的密度分布圖Fig.9 Density distribution of samples in the virtual screening space

如圖9 所示可以看到，篩選空間中的分子的預(yù)測(cè)晶體密度全部高于1.65 g·cm-3，且預(yù)測(cè)晶體密處于1.8～1.9 g·cm-3范圍的分子占到了大多數(shù)，晶體密度接近2.0 g·cm-3的分子數(shù)量達(dá)到了28。這一結(jié)果表明，采用預(yù)篩選高密度骨架的方式構(gòu)建的含能分子篩選空間，其分子樣本具有顯著的高晶體密度特征。這一結(jié)果驗(yàn)證了圖5 中的分析結(jié)果，即分子骨架等效密度與晶體密度具有正相關(guān)性，預(yù)篩選高密度分子骨架可以有效提升高密度含能分子的篩選成功率。

為了評(píng)估本研究設(shè)計(jì)的含能分子的綜合性能，選取密度排名前100 的分子，由EM-Studio 計(jì)算其生成焓、爆速、爆壓和BDE。為了直觀體現(xiàn)被篩選分子的綜合性能，采用常見的含能分子TNT、RDX 和TATB 作為指標(biāo)分子，對(duì)各種性質(zhì)進(jìn)行了對(duì)比分析（圖10）。

圖10 候選含能分子的性能，（a）爆壓，（b）爆速，（c）BDE 與（d）晶體密度（紅色線為TATB 的對(duì)應(yīng)性質(zhì)，綠色線為TNT 的對(duì)應(yīng)性質(zhì)，藍(lán)色線為RDX 的對(duì)應(yīng)性質(zhì)）。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 中我們可以看到，有超過半數(shù)的分子爆壓超過TATB（圖a），并且爆速整體超過TATB，有部分分子爆速甚至超過RDX（圖b），在BDE 上，大部分分子超過RDX，過半分子超過TNT（圖c），全部分子密度超過TNT，并且大部分超過RDX（圖d），這說明通過密度篩選出來的分子的爆轟能力普遍大于TATB 和TNT，并且個(gè)別分子甚至超過了RDX，可能出現(xiàn)的原因是因?yàn)閷?shí)驗(yàn)篩選出的分子的氮含量都比較高，且環(huán)的數(shù)量都大于等于2，這類分子被認(rèn)為［49］環(huán)內(nèi)的N—N 鍵和C—N 鍵在爆炸時(shí)會(huì)釋放出大量的能量，可有效提升爆轟性能。從圖10c 中可以看到有接近一半的分子穩(wěn)定性優(yōu)于RDX，其原因可能是因?yàn)榉肿庸羌苤写罅亢羞辉?、四嗪、四唑和三唑等富氮雜環(huán)，研究表明［50］，這些分子具有良好的熱穩(wěn)定性，已經(jīng)成為部分研究設(shè)計(jì)含能分子的有效結(jié)構(gòu)單元。

在篩選分子的時(shí)候?qū)⒈僖訰DX 為基準(zhǔn)，BDE 以TNT 為基準(zhǔn)進(jìn)行對(duì)比，以期得到爆轟性能強(qiáng)于RDX，且穩(wěn)定性強(qiáng)于TNT 的含能分子（圖11）。

如圖11 所示，在BDE 與爆速同時(shí)較高的分子數(shù)量較少，爆速超過RDX，BDE 超過TNT 的分子僅僅只有6 個(gè)，進(jìn)一步驗(yàn)證了含能分子的能量與安全性之間可能存在一些矛盾，高能量的分子往往穩(wěn)定性較差，在設(shè)計(jì)分子的時(shí)候應(yīng)高要綜合考慮爆轟性能與穩(wěn)定性的關(guān)系。研究表明［3］，在分子水平上，能量和穩(wěn)定性的矛盾尤為突出，但在晶體和混合物層次，能量和安全性的矛盾可以得到很大程度上的緩解。因此隨著后期的晶體工程和復(fù)合工藝的發(fā)展，可能有希望得到越來越多的高能量高穩(wěn)定性的含能分子。這6個(gè)分子的環(huán)中的化學(xué)鍵均為大π 鍵，并且多個(gè)環(huán)之間僅有一條化學(xué)鍵相連，或許這樣的結(jié)構(gòu)能在保證分子化學(xué)穩(wěn)定性的基礎(chǔ)之上，較大的提升分子在爆轟時(shí)釋放出來的能量，后續(xù)在設(shè)計(jì)分子和選擇骨架的時(shí)候，應(yīng)多選擇此類骨架和分子。具體分子性質(zhì)如表2 所示。

圖11 Ⅰ區(qū)為候選分子爆速與BDE 的權(quán)衡，Ⅱ區(qū)為在Ⅰ區(qū)中右上角爆速大于RDX 且BDE 大于TNT 的分子結(jié)構(gòu)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 中可以看到生成的分子在爆轟性質(zhì)上基本超過TATB，在部分性質(zhì)上甚至超過RDX，化學(xué)穩(wěn)定性均強(qiáng)于RDX，因此可以嘗試將本次實(shí)驗(yàn)設(shè)計(jì)的化合物進(jìn)行進(jìn)一步的研究。值得關(guān)注的是，本研究設(shè)計(jì)的分子在穩(wěn)定性上仍無法超過TATB。TATB 在結(jié)構(gòu)上具有獨(dú)有的特點(diǎn)，其3 個(gè)硝基和3 個(gè)氨基是相間排列在苯環(huán)上的，使整個(gè)分子體系的電荷分布高度對(duì)稱，相間分布的氨基和硝基通過強(qiáng)的分子內(nèi)氫鍵和靜電作用在分子外圍形成非常穩(wěn)定的環(huán)狀結(jié)構(gòu)［53］。從當(dāng)前的研究來看，設(shè)計(jì)出類似于TATB 這樣完美的分子結(jié)構(gòu)非常困難，這可能是探索新型高穩(wěn)定性含能分子需要面對(duì)的挑戰(zhàn)。

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

3 結(jié)論

（1）含能分子的晶體密度與骨架密度成正相關(guān)，基于高密度骨架可構(gòu)建具有潛在高密度的待篩選分子。

（2）通過分子骨架預(yù)篩選獲得200 個(gè)高密度的稠環(huán)分子骨架，用硝基和氨基作為取代基，可以組合出7000 多個(gè)潛在的含能分子，且均具備高密度的潛質(zhì)。

（3）基于CCDC 數(shù)據(jù)庫中的硝基化合物為數(shù)據(jù)集，通過圖神經(jīng)網(wǎng)絡(luò)訓(xùn)練出晶體密度預(yù)測(cè)模型，其決定系數(shù)R2達(dá)到0.924，且泛化性良好。

（4）基于深度學(xué)習(xí)、量子化學(xué)計(jì)算和爆轟產(chǎn)物狀態(tài)方程等方法預(yù)測(cè)晶體密度、生成焓、爆轟性能和化學(xué)穩(wěn)定性，從而由性能排序篩選出能量水平優(yōu)于RDX，穩(wěn)定性優(yōu)于TNT 的新型含能分子6 個(gè)。

（5）通過預(yù)篩選分子骨架可以有效提升虛擬篩選空間的總體性能，在此基礎(chǔ)上可借助高通量計(jì)算與深度學(xué)習(xí)實(shí)現(xiàn)含能分子的高效設(shè)計(jì)。