机器学习势在含能材料分子模拟中的研究进展

常晓雅,文明杰,张 迪,王永锦,初庆钊,朱 通,陈东平

(1.北京理工大学 爆炸科学与技术重点实验室,北京 100081;2.上海市分子治疗与新药开发工程研究中心, 华东师范大学 化学与分子工程学院,上海 200062)

引 言

含能材料通常指在受到外界刺激时能快速分解并释放巨大能量的化合物或混合物[1-2],被广泛应用于军事燃烧及爆破、火箭导弹发射、推进剂等领域。在推进剂领域,含能材料以燃烧的形式释放能量,实现了从化学能向动能的有效转化,成为先进导弹武器发动机的动力来源,是推动武器装备系统发展的支撑和制约技术之一[3-4];在火炸药领域,含能材料在极短的时间内释放大量的热能和气体,以爆轰形式作功形成冲击波,是武器装备完成发射和高效毁伤的能源材料。作为现代武器系统的主要能量来源,含能材料燃烧与爆炸性能直接影响武器装备的打击范围和毁伤能力,是先进武器装备优化设计的关键发展方向。

含能材料释放的能量主要来源于原子核外层电子在反应过程中转移所释放的能量[5],其宏观特性本质上是由其微观分子结构所决定的。例如,含能材料爆轰反应放热量取决于产物分子的生成热;含能材料的感度与分子结构、分子间相互作用以及反应活化能等性质有关;各组分之间的相容性与分子间的相互作用有关[6-7]。因此,在微观结构层面上研究含能分子结构燃烧特性和力学稳定性对于新型含能材料的设计至关重要。分子模拟作为原子尺度的研究方法,可以用来预测含能材料的力学行为和化学反应过程,揭示分子结构与宏观燃烧性能之间的关系,从而为含能材料的理论设计提供科学依据。在分子模拟计算中,描述微观体系中不同分子构型与能量的映射关系被称为势能面。分子模拟的实质是在被模拟材料分子体系的势能面上进行动力学研究,例如搜寻分子的最低能量构形或者计算化学反应速率。其中构建势能面的计算方法主要包括第一性原理计算和分子动力学。第一性原理计算是通过求解薛定谔方程来描述粒子的电子结构,在不需要任何先验参数情况下,准确计算分子、原子的理化性质,计算精度较高。在空间尺度较小且精度要求高的场景,比如少量分子体系的化学反应路径研究、动力学参数计算等,第一性原理计算表现出独特的优势。然而高精度的计算水平伴随着高昂的计算成本,时间复杂度往往是O(N3)甚至更高,现有的计算能力无法承载超过1 000个原子的分子体系。

分子动力学不考虑电子运动的自由度,电子效应体现在原子间相互作用的势函数中,原子的运动遵循牛顿第二定律,并据此实现目标体系在给定初始条件下随时间的动态演化。分子模拟的精度取决于势函数对真实势能面重建的精度。势函数的数学形式和分子各部分对能量贡献的可拓展性等均对势能面重建精度有影响[8]。非反应力场假设的函数形式简单,如施加两体或三体截断,对于材料的力学性质描述较为准确,其研究尺度可达上百万甚至千万原子的体系[9]。但这样的经验力场并不适用于研究电子转移、原子变价的化学反应过程,无法用于模拟含能材料燃烧与爆炸等复杂反应过程。反应力场是基于键级来描述体系的键合作用,键级随着化学反应的发生自发进行调整[10],包括Airebo[11]、COMB3[12]和ReaxFF等力场[13-14]。其中,应用最广泛的是Van Duin等[13-14]提出的ReaxFF力场,通过对键能、范德华能等多个能量项进行参数拟合来构建势能面,目前在含能材料点火、燃烧、冲击等方面均有大量研究[15-20]。利用反应力场进行分子模拟,结合模拟轨迹分析,可得到反应中间体、反应产物等随时间的变化,从而提取其主要的基元反应、构建详细反应机理。精确的力场模型是探究含能材料反应机理的基础。此外,并不是所有含能材料体系均已开发ReaxFF力场,以AP为例,尽管已经广泛验证其作为氧化剂可以增加含能材料的能量释放和反应活性,但是并未有前人开发相应的ReaxFF力场,因此目前并没有关于AP燃烧过程的分子模拟研究。

第一性原理计算精度高,但计算成本大;分子动力学模拟效率较高,但是精度较低。两种方法均无法兼顾精度与效率。如何高效率、高精度地重建势能面是研究人员面临已久的难题。近年来,基于机器学习的势函数模型有望解决传统计算方法的“卡脖子”问题。该方法绕过理论模型精确解释的难题,通过完全数值化对抽象的物理规律和复杂的化学过程之间的潜在关系进行推断,从而准确高效地构建反应体系的势能面,实现多种复合含能体系的理论计算,建立分子结构与燃烧性能的联系。

本文归纳了机器学习势函数模型搭建方案及训练集搭建策略,通过对高精度的量子化学势能面进行近似,成功应用于含能材料燃烧爆炸方面的研究,包括常用的RDX、高氮炸药CL-20、新型硝胺炸药ICM-102、强氧化剂AP、高能颗粒(铝和硼)等物质,并综述了机器学习势在碳氢燃料燃烧中的研究进展。最后总结了机器学习精度和效率方面的优势,并对其发展趋势提出展望。

1 机器学习势

1.1 机器学习势模型

精确构造含能材料的势能面是开发机器学习势的基本要求,也是研究其微观燃烧机制的基础。机器学习势从第一性原理数据出发,用高维函数来描述分子结构与能量之间的关系,进而拟合势能面,以较低的计算成本实现对目标体系的分子模拟计算。一个合格的机器学习势模型需要至少满足广延性、对称不变性和光滑性[9]。体系的广延性是指体系总能量和原子数成正比,故机器学习势的能量分解形式为:

(1)

(2)

为了达到上述要求,通常引入描述符D(Ri)的概念对原子三维坐标进行映射。对于进行等价对称操作的原子坐标,其描述符相同,即:

D(Ri)=D(URi),∀U∈u

(3)

从而将原子能量贡献改写为:

(4)

其中F为拟合方法。将原子坐标映射得到的结构描述符作为输入,通过不同的拟合方法得到分子结构和能量的对应关系,即机器学习势模型。模型的光滑性由拟合方法的光滑性和描述符共同保证。本节简要介绍机器学习势的拟合方法和描述符构造方法。

目前主流的机器学习拟合方法包括核方法和神经网络方法。其中,核方法通过测量描述符与数据集中描述符的相似程度来获得系统的能量或力。2010年,Csanyi等[21]基于核方法提出高斯近似势(Gaussian Approximation Potential, GAP)模型,通过对原子化学环境相似性进行度量,采用高斯回归过程对势能面进行内插,从而得到原子能量、受力的最优估值[22]。分子构型越相似,在势能面上的位置就越接近,构型的势能也应当越相近[23]。这种方法并未预设函数形式,训练过程人工干预较少,且可解释性较强。目前在材料领域已经发表了大量的相关研究工作,如钙钛矿[24]、过渡态金属[25-26]、硅烯[27]、半导体[28]、合金[29]等。2018年以来,基于核方法的(对称)梯度域机器学习((Symmetric)Gradient Domain Machine Learning, GDML)[30]方法快速发展。与高斯近似势不同,GDML模型在初始数据集中仅加入第一性原理计算的原子受力,通过对原子受力进行积分得到体系总能量[31]。引入受力项的模型比基于能量的模型自由度更高、参数更多、形式更灵活。

神经网络势函数(Neural Network Potential, NNP)模型是用一个多参数的非线性函数来拟合势能面。早在1995年,Blank等[32]利用前馈神经网络拟合势能面,研究CO在Ni表面以及H2在Si表面的扩散问题。该模型的输入坐标并不满足平移、旋转和交换对称不变性,故前馈神经网络仅仅应用于低维体系。2007年,Behler和Parrinello等[33]提出了高维神经网络(High-Dimensional Neural Network, HDNN)模型, 采用以原子为中心的对称函数矢量来描述原子的邻域环境,每个原子用一个单独的前馈神经网络表示,从而预测每个原子对体系全局能量的贡献。这样的神经网络预测能力更强,但需要反复调节和拟合大量的网络参数, 且需要充足的训练样本。一般来说,神经网络模型的深度比其宽度更受青睐,即隐藏层的层数比每层的节点数更重要[34]。这是因为更深的神经网络具有更好的非线性表达能力,可学习更复杂的变化,从而拟合更复杂的特征。换言之,包含多个连续隐藏层的深度神经网络的高维非线性函数拟合能力更为强大。2018年,鄂维南、张林峰等[35-36]基于深度神经网络提出了深度势(Deep Potential, DP)模型。同年,该团队开发了开源工具DeePMD-kit[37]来支持DP模型的开发工作。目前DP模型已被用于材料、力学、化学、生物等各个领域的研究,如电池[38-39]、催化[40-41]、蛋白质[42-43]、含能材料[44-46]等。以DP模型为代表的机器学习势模型,通过机器学习算法加速和大规模并行优化,可在保持第一性原理计算精度的同时实现上亿原子的分子动力学模拟,将精确的物理建模带入了更大尺度的材料分子模拟[47]。

机器学习势通过构造结构和能量之间的映射关系来重建势能面,因此采用何种描述符来提取分子结构特征进而用于训练机器学习模型对于力场开发的精度和效率至关重要。Behler和Parrinello等[33]构造的原子中心对称函数(Atom-Centered Symmetry Function, ACSF)满足物理对称性以及相对于移动原子的可微性,被广泛集成于多个软件中[48-49]。该方法以原子为中心,用径向函数和角向函数将邻域原子的化学信息编码为复杂的“多体”形式,以完整表示原子的局域环境,并引入一个平滑函数,使其在截断半径处平滑过渡到零,以保证原子运动的平滑性和连续性。值得注意的是,ACSF方法给出的是构型的绝对描述,非常适用于在高维神经网络势中进行原子化学环境描述。

Bartok等[50]并未试图建立局域构型的描述,而是直接建立两个局域构型之间距离的对称性度量,即平滑重叠原子位置(Smooth Overlap of Atomic Position, SOAP)。这种方法的主要依据是局域原子密度,将两个局域构型之间的“距离”定义为两个局域密度之间的交叠程度,并将核原子与其余原子的相对位置表示为径向项和角向项的乘积。对于原子数量较多的团簇,SOAP描述符具有独特的优势,通常能得到更准确和稳定的势能面。然而,SOAP给出的并非是某个构型的绝对描述,而是构型之间的距离,因此并不适合用于深度神经网络,而非常适用于核方法[9]。同样地,库仑矩阵(Coulomb matrix)[51]描述并非针对局域原子环境,而是针对体系中所有的原子,将体系中N个原子的三维坐标转换为一个N×N的矩阵,多用于开放边界的体系,而非周期性体系的描述。库伦矩阵保证了坐标的平移和旋转不变性,对库伦矩阵排序后可以保证交换对称性。

描述符构造是开发机器学习势的基础,只有将描述符和拟合方法有效结合,才能最大程度地发挥机器学习势的优势。需要注意的是,关于何种分子描述符在机器学习势模型中表现最优、从原子坐标中提取哪种特征,目前学界并无定论[8]。

1.2 数据集搭建方法

1.2.1 采样方法

机器学习势的开发过程中需要首先利用高精度的量子化学计算方法来构建体系的势能面。从效率和精度的角度出发,通常采用领域中广泛使用的密度泛函(DFT)方法。由此可见机器学习势的精度上限是量子化学数据集的精度。考虑到DFT的计算成本,完全覆盖复杂体系高维全局势能面也是不切实际的,因此开发机器学习势的核心是搭建合理的数据集。数据集的充分性、多样性对于机器学习势的精度至关重要。训练集需要尽可能地覆盖目标体系的势能面,对不同含能材料组分、多种原子局域环境进行批量采样,以得到精度高、拓展性强、普适性高的分子动力学机器学习势函数模型。常用的势能面采样方法(见图1)包括结构微扰法、分子动力学(MD)模拟方法和增强采样方法。其中结构微扰法是对原子坐标施加结构缩放、错位和随机扰动,从而得到大量随机分子结构;MD采样方法采用现有的反应力场、经验力场等进行分子动力学计算,得到分子的运动轨迹和反应路径;增强采样法则是通过人为施加约束,降低化学反应能垒,提高反应发生的概率。

图1 常见的采样方法:结构微扰法、MD采样法、增强采样法和虚拟现实MD采样法[45]

含能材料组分复杂,通常涉及到3种以上的化学组分,5种以上化学元素;且多应用于高温(>3 000 K)、高压(>10 GPa)等极端场景,其反应过程涉及到气、液、固3个相态及两相甚至多相之间的非均相反应,这些条件决定了含能材料机器学习势的开发难度远大于其他材料,往往需要更大规模的数据集来覆盖可能存在的热力学条件。因此在进行势能面采样时,需要有针对性地构造数据集。例如,在高温条件下,分子化学键断裂生成大量自由基,原子局域环境的复杂度显著增加,因此初始数据集中需包括已知的反应路径、过渡态结构和反应产物。为了描述含能材料的高压构象空间,需要在数据集中引入不同压力状态的结构构象,即材料的状态方程(压力—体积曲线),以描述含能材料在极端场景的力学特性。

事实上,在传统分子动力学模拟的时间尺度内,体系很难跨越较高的反应能垒,在相空间内的采样效率具有明显的局限性,使用增强采样方法能够更高效地对稀有事件(化学反应)实现采样。然而,经典的增强采样方法需要预先选定反应坐标,通过施加偏置势或偏置力使得体系沿着既定的方向发生反应,考虑到含能材料复杂的分子结构和苛刻的热力学条件,潜在反应坐标数量指数性增长,因此经典的增强采样方法会存在反应坐标无法合理确定的困境,以至于对潜在分解反应路径、过渡态结构等采样不足。近年,本课题组[52]开发了一款虚拟现实(virtual reality, VR)分子动力学模拟软件MANTA,为含能材料数据集的搭建提供了新的采样工具。研究人员可以借助VR技术并凭借化学直觉对目标反应路径开展实时交互,快速地对特定的反应过程势能面进行采样,将反应构型实时加入到训练集中,是一种极具应用潜力的辅助采样工具。总之,势能面采样是开发机器学习势的关键步骤,在训练过程中应结合上述手段,尽可能完备地描述化学反应的势能面空间,以提高势函数的准确性。

1.2.2 主动学习

主动学习是指在机器学习势训练过程中,对目标体系的相空间进行不断地探索和补充,实现自动拓展数据集的一种训练方法[53]。对于复杂体系,主动强化学习是一种十分合适的数据挖掘、构型探索、势能面拓展方法。以DP-GEN[54]为例,这一工具是在DeePMD-kit的基础上开发的同步学习工作流,其工作流程如图2所示。

图2 DP-GEN工作流程

针对初始数据集,调用DeePMD-kit工具,利用不同的神经网络初始化设置,训练多个机器学习势函数模型。随机选取其中一个势函数模型进行分子动力学模拟(调用LAMMPS等[55]软件),对运动轨迹中的构型能量和受力进行评估。若该模型预测构型的能量或受力与其他几个机器学习势函数模型有明显偏差,则表示当前模型对于该构型的预测精度较差,需要补充到初始数据集中从而进一步优化势函数模型。DP-GEN引入模型偏差来判断多个机器学习势间的误差,通过设定误差上下限来挑选候选构型。随机挑选其中部分构型进行DFT计算得到新的数据,并加入到现有的数据集中重新训练,从而得到精度比之前更高的机器学习势模型。上述流程对数据集进行补充挖掘并得到新的机器学习势模型,被称为完成了一轮迭代。

采用DP-GEN方案进行多次迭代,每轮迭代可以改变系综设置,完成在目标温度、压力或体积条件下的批量采样和定量筛选,实现主动学习,旨在有效覆盖所需样本空间。高效的迭代训练过程设计能让势能面覆盖区域从局域空间平滑稳定地拓展到整个构象空间,例如从短的模拟时间逐步增加为长时间模拟,从低温模拟逐步扩展为高温模拟,从低压模拟逐步增加为高压模拟,从非反应体系过渡至反应体系。

2 基于机器学习势的含能材料分子模拟

依据材料功能,含能材料组分一般可分为五大类,如图3所示。

图3 含能材料常见组分

其中,硝胺类是指以碳、氢、氧、氮等元素构建的具有爆炸性基团的化合物,是含能材料的主体部分[5]。硝胺类组分以环状和笼状的硝胺硝基化合物为主,如环三亚甲基三硝胺(RDX)、环四亚甲基四硝胺(HMX)、六硝基六氮杂异伍兹烷(CL-20)等。这类分子受限于反应产物的生成热,导致其能量密度偏低。高能颗粒(铝、硼、镁等)凭借较高的燃烧焓被广泛用做还原剂以提高含能材料配方体系的能量密度,调节其能量输出特性,是强有力的动力来源和毁伤手段[56-58]。然而,高能颗粒的氧化层往往会限制其能量释放速率,导致较长的点火延迟时间。在含能复合材料中加入聚偏二氟乙烯(PVDF)、高氯酸铵(AP)等强氧化剂有助于完全发挥纳米颗粒的增益作用,其燃烧和热分解性能极大改变了传统含能材料的释能机制[59-60]。为了降低含能材料感度,提高其力学性能,组分配方中通常包含少量黏结剂作为辅料。黏结剂组分多为高分子材料,如端羟基聚丁二烯(HTPB)、聚叠氮缩水甘油醚(GAP),这些高分子材料同样包含一定含能基团以提高整体配方的能量密度。此外,近年来一部分高能量特性、热稳定性、机械感度及力学性能优异的催化剂在含能材料领域中得到了广泛的应用,如氧化石墨烯(GO)、铜(Cu)等[61-62]。

本章将介绍机器学习势在含能材料分子模拟方面的研究进展,涉及到硝胺类化学物(RDX/CL-20/ICM-102)、氧化剂(AP)、金属添加剂(Al/B)等组分的反应动力学研究。同时,包含了机器学习势在碳氢燃料燃烧分子模拟中的进展。上述研究所采用的机器学习模型均为Zhang等[37]提出的DP模型,简记为NNP。

2.1 硝胺类含能材料

2.1.1 RDX

黑索金(RDX)是一种典型的含能材料,稳定性好、爆炸威力大,被广泛应用于火炸药领域。Chu等[44]通过CP2K计算[63]将PBE/DZVP级别的第一性原理数据作为训练集,基于DP策略开发机器学习势模型对RDX热解的反应动力学进行研究,并与ReaxFF[16]的结果进行对比。两种方法得到的不同温度下RDX分子数量变化趋势相同,即体系分解速率均随温度升高而增加;但变化快慢却大不相同,NNP模型预测的分解速率与第一性原理计算的结果一致,约为ReaxFF力场结果的两倍。这表明开发的机器学习势函数能准确描述RDX的热解反应,而ReaxFF反应力场低估了分解反应的反应速率。此外,NNP模型和ReaxFF力场通过阿伦尼乌斯公式得到的RDX热解反应活化能分别为105.6 kJ/mol和108.8 kJ/mol,与文献值119.7 kJ/mol接近,即两种力场模型均能较为准确地预测RDX分子发生单分子分解反应的活化能。

RDX热分解的中间体和平衡产物主要包括NO2、NO、HNO2、N2、H2O、CO2和CO等。机器学习势函数以第一性原理的精度成功复现文献结果,即RDX受热时N—NO2键会首先发生均相断裂而生成NO2,中间产物NO2进一步发生反应生成NO和H2O分子。当达到平衡浓度时,体系中N2含量最多,其次是H2O和CO2。与ReaxFF模型的结果相比,机器学习势预测到更多的NO分子生成量、更快的H2O生成率以及更高的N2平衡浓度,与文献结果一致[64]。通过对分子动力学模拟的轨迹进行分析可以得到RDX热解的详细反应路径,如图4所示。RDX的分解机制可分为两个阶段:(1)N—N键断裂生成NO2分子和RDR中间体;(2)中间体RDR进一步发生N—N键断裂、开环反应、或氢加成反应。其中,关键中间体NO2发生氢加成反应会生成HONO,进而生成NO、N2、OH和H2O等物质。

图4 RDX分解的中间产物及详细反应路径[44]

近邻分子间的相互作用在很大程度上会影响甚至改变含能材料的分解路径[6]。不同密度的RDX分子和RDR中间体发生N—N均裂、氢提取和开环反应的通量并不相同。在高密度情况下,由于空间位阻效应,近邻分子会阻碍RDX分子及其中间体的开环反应;随着密度降低,近邻分子的影响减弱,开环反应通量增加,氢提取反应通量减少,RDX的热解机制与孤立分子的第一性原理结果相近。上述研究结果在保持第一性原理计算精度的同时,从原子尺度阐明了RDX的反应路径,建立了RDX热分解详细反应动力学机制。

2.1.2 CL-20

第三代高爆炸药CL-20是能量和密度最高的单质炸药之一,在氧平衡、爆速、爆压和堆积密度等方面均优于其他含能材料。然而,对机械刺激的高敏感性限制了其应用前景。近年来,共晶技术在含能材料领域已成为一种广泛应用的策略,以更好地平衡感度和能量输出之间的关系。现有的实验和理论工作表明,CL-20与其他含能材料共晶结合可以提高冲击稳定性和热稳定性。尽管已经有许多关于CL-20的共晶材料被设计出来,但提高稳定性的具体微观机理仍有待进一步揭示。

Cao等[46]构建了β-CL-20和CL-20/TNT共晶体系PBE/DZVP级别的神经网络势能函数,并在此基础上开展分子模拟研究以揭示CL-20/TNT共晶稳定性增强背后的分子机制。为了深入探索β-CL-20和CL-20/TNT共晶的热解过程,基于机器学习势进行了分子动力学升温模拟。对于β-CL-20体系,CL-20分子在温度升高到大约1 200 K时开始分解;而在CL-20/TNT共晶体系中,CL-20分子发生分解的温度高达1 400 K左右。通过比较起爆温度,可以发现共晶CL-20/TNT的热稳定性明显提高。在不同的升温速率下,也能观察到相似的结果,即在CL-20/TNT共晶体系中的起爆温度总比在β-CL-20体系中的要高。这是由于共晶体系会稳定生成分子间氢键,而在单质体系中并没有氢键的出现。分子间的氢键相互作用可以有效延迟CL-20分子的裂解,降低活性物种的反应速率,这一结果与张朝阳等[65]对含能材料晶体结构的研究结果一致。通过对分子动力学模拟轨迹进行分析,得到单质CL-20分子主要的热分解路径,如图5所示。与之前的研究一致,CL-20分子热解以N—NO2键的断裂为主,脱去NO2后的中间体进一步发生N—N断键和C—C断键,生成产物NO2和CxNy小分子。在共晶体系中,也观察到了类似的初始热解路径,但是速率较低。TNT中的自由基与游离的NO2分子会发生进一步反应,如TNT中的氢自由基与游离的NO2分子结合生成HONO,TNT中的甲基自由基与NO2结合生成甲醛(HCHO)和NO。由于TNT中的C—N键和六元环相对稳定,因此在模拟过程中很难发生NO2的分离和开环反应。值得注意的是,从CL-20中分离出的NO2与TNT的反应进一步减缓了初始阶段NO2的生成速度,阻碍了CL-20分子的进一步分解。因此,TNT可以被视为缓冲碎片,以降低活性中间体与CL-20分子碰撞的概率。上述工作从原子尺度揭示了CL-20/TNT共晶体稳定性增强背后的物理化学机制。

图5 CL-20热解的初始反应路径[46]

2.1.3 ICM-102

新型含能材料的设计是提升我国先进武器装备综合性能的有力支撑。张庆华等[66]将“材料基因组工程”的研发模式引入到了含能材料领域,针对含能材料能量与安全匹配性差的问题,运用计算机辅助分子高通量设计与筛选的手段,成功研发了一种高能钝感炸药分子ICM-102(C4H6N6O4)。其能量水平高于低感炸药LLM-105,感度与钝感炸药TATB相当[67],实现了炸药能量与安全性的平衡。Chu等[45]基于虚拟现实增强采样(VRMD)方案构建了PBE/DZVP级别的数据集,并开发机器学习力场来描述ICM-102的反应动力学。其中,第一性原理计算采样的构型与初始构型接近,随温度的升高,构型在势能面上的分布区域扩大,但主要集中于低能区域,导致反应构型数量有限,不能完整描述整个反应过程。而VRMD方案通过手柄施加外力,并依靠化学直觉对经验的反应位点、反应路径进行实时交互采样(见图2)来构建数据集,可以直接对断键反应等稀有事件进行直接采样,有效提高了高能垒反应过程的采样效率。

数据集的质量对于机器学习模型的性能至关重要。通过主成分分析可提取采样结构的主要特征并进行数据降维,从而将采样到的结构划分为晶体构型、部分分解构型和气态产物,如图6所示。沿着总反应发生的方向,即主成分PC1增加的方向,两种采样方案得到的构型质量均逐渐减少,但是对应的数密度却明显不同。由于反应能垒的存在,直接分子动力学采样的数据集范围明显集中在低能势阱,在高能区域的采样较少,即得到的构型大多并未涉及到反应的发生。而VRMD采样方案可针对能垒较高的反应过程进行采样,作为AIMD采样方式的补充,极大地拓展了势能面的采样空间。此外,单一的AIMD采样方案得到的模型能量精度为0.014 eV/atom,两种方法相结合得到的模型偏差仅有0.006 eV/atom。这表明VRMD采样方案可以作为补充手段,来构造更为均衡、更丰富的训练集,从而得到更精准、更有效的势函数模型。

图6 AIMD和VRMD采样方案的主成分分析[45]

通过高温分子模拟,Chu等[45]发现ICM-102热分解的主要中间体和产物包括C6H5N6O4、H2O、NO、N2、CHNO和CO2。利用聚类分析的方法可以将所有出现过的分子构型分为八类:前两类分子表示ICM-102及其失去一个小分子基团后的产物;第三、四类分子表示ICM-102分子经过开环反应的产物;第五、六、七类分子表示开环反应后进一步生成的含碳或含氮的分子;第八类则是反应过程中的水、一氧化氮和氮气等小分子气体。以热解反应网络为基础,可直接依据通量分析推演任意中间物质的详细反应网络(见图7)。

2.2 AP

AP的有效氧含量高、机械感度低、相容性好,爆炸过程时产生的气体多、无固体残渣、成本低,是固体推进剂和火炸药等最重要的强氧化剂[68-69]。AP的热分解是一个复杂的物理化学过程,涉及到吸热相变和放热反应。目前AP的热分解研究以传统实验为主,如使用热重分析(TGA)、差示扫描量热法(DSC)和质谱分析(MS)等来确定分解速率和中间产物。传统的实验手段并不能在极小的时间、空间尺度上直接观测到AP复杂的反应机制。而AP反应力场的缺失制约了研究者对其热分解机理的微观研究,是火炸药领域中的“卡脖子问题”。Chu等[70]基于第一性原理计算(PBE/DZVP级别)首次建立了AP的机器学习势,对其热分解过程进行分子动力学模拟,从原子尺度成功提取AP分解的动力学演化过程,填补了目前含能材料燃烧机理缺失的重要部分。通过分子模拟复现TGA/DSC实验,将AP晶体从常温加热至3 000 K,所开发的机器学习势可描述AP热分解过程中的质量变化和热流速率,即从正交相到立方相的吸热相变,以及随后的放热分解反应。

AP分子由NH4和ClO4基团构成,通过直接分子动力学模拟,机器学习势模型预测其分解过程中生成的中间体和反应产物多达三万种。其中,反应的中间体主要是NH3、ClO2、ClO3、HClO4。随着反应的进行,最后的平衡产物以H2O、NO、O2、Cl2、N2和HCl为主。通过提取热分解过程中的主要物质和生成路径得到AP反应网络,如图8所示。

图8 AP分解的产物数量变化和主要反应路径[70]

考虑到AP中的强氧化性元素较多,据此可以将整个分解动力学划分为AP的初始分解、氢化学、氮化学和氯化学。初始分解过程以质子转移为主,NH4中的氢自由基会转移到HClO4中并生成NH3,或与ClO4中的O自由基结合生成ClO3和OH。在初级分解之后,氮化学、氢化学和氯化学反应同时进行。特别指出的是,上述分子模拟预测到一些实验未曾报道过的双分子反应,例如NH4+NH2→2NH3和ClO4+NH3→HClO4+NH2。这些气相物种和AP之间的双分子反应可能会影响AP的分解过程。上述研究以第一性原理计算的精度首次揭示了AP的详细反应机制,克服了AP反应力场缺失对基础研究的制约,为建立复合含能体系的反应动力学模型提供了新的机遇。

2.3 高能颗粒(Al、B)

金属颗粒能量密度大、燃烧热值高,是极具吸引力的含能材料组分,能显著提高固冲发动机的性能,增大导弹武器的射程,提高战术武器的作战效能和毁伤效益[40, 71-72]。金属作为燃料和燃料添加剂可以减少燃烧不稳定性、提高燃烧速率、降低爆震敏感性[57, 73]。其中应用最广的金属颗粒是铝,其燃烧热值高、能量密度大且储量大,是重要的含能材料添加剂。近期本课题组基于PBE/DZVP级别的第一性原理数据建立机器学习势模型,训练集中包括Al、Al2O3、界面以及氧气扩散等数据。所开发的模型可以复现晶体铝的状态方程曲线(见图9),计算结果与DFT预测基本一致。模型在极端的压缩和拉伸情况下仍然保持着良好的性能表现,因此可以用于研究铝在冲击波下的力学响应行为。值得注意的是,尽管训练集中仅包含了单位原子体积在10～25 Å3/atom的数据(图中虚线区域的数据),机器学习势模型可以准确地预测7～32 Å3/atom范围内的状态方程,这证明了机器学习势函数具有一定的外推性能。

图9 铝的状态方程曲线

在所有金属和类金属中,除有毒性的铍外,硼的体积热值和质量热值最高,分别为138 kJ/cm3和59 kJ/g[74-75]。尽管硼的理论燃烧性能优异,但是硼的熔点(2 348 K)和沸点(4 043 K)都很高,导致颗粒难以点火,限制了硼粉的实际应用。此前Liu等[76]、Wang等[77]采用分子模拟的方式研究了硼颗粒的微观燃烧机理,然而这些工作使用的势函数是基于氨硼烷[78]开发的,并未针对硼燃烧进行优化。这里我们提出一个问题:一个力场模型的开发数据集中并不包括硼单质的数据,这样的模型可直接用于硼燃烧的模型吗?势函数模型是针对目标体系特定的应用场景而开发的。一般而言,训练集只是全局势能面上的局部空间,在势能面上的外推性能非常局限。例如,同样含有C/H/O/N元素的反应力场可以描述RDX的热解反应,却不能直接用于CL-20的模拟,会导致预测误差迅速升高。因此,在进行分子动力学模拟之前,需要对势函数开发的训练集进行分析,确保数据集覆盖了关注的势能面区域,避免超出训练集的外推空间,并进行物理化学参数验证,如晶格参数、状态方程、融点、比热容等。

Chang等[79]利用CP2K构建了PBE/DZVP级别的量子化学数据集,开发了用于描述硼和硼氧化合物的DP模型。开发的DP模型对能量、受力、状态方程、晶格参数和径向分布函数等进行验证,机器学习力场的精度均和第一性原理计算一致。相比之下,针对氨硼烷燃烧开发的ReaxFF力场预测结果非常差[78]。例如,文献中晶体硼的单位原子体积为7.25 Å3/atom,机器学习预测值为7.01 Å3/atom,而ReaxFF力场的预测值高达12.58 Å3/atom,显著高估了晶体硼的平衡结构。以上这种明显误差主要来源于ReaxFF力场的滥用,考虑到上述工作的ReaxFF是针对氨硼烷燃烧体系,在其开发过程中不含任何硼单质晶体结构,因此这套ReaxFF力场参数无法支持硼燃烧的分子动力学模拟。另外,这直接表明研究人员应当避免任意拓展ReaxFF力场的应用范围。

Chang等[79]采用DP模型模拟了纳米硼颗粒的融化行为。通过计算颗粒全局Lindemann指数和势能变化,得到不同粒径硼颗粒的熔点,归一化后的熔点如图10所示。对于3 nm硼颗粒,其融点高达1 992 K,而同样粒径的铝镁颗粒在该温度下已完全融化并蒸发到气相环境中[80-81]。这表明纳米硼颗粒的燃烧机制仍然是以表面反应为主。通过分析原子Lindemann系数演化,作者发现硼融化的过程包括初始加热、表面预融、内核预融、完全融化和继续升温等五个阶段。该融化过程遵循Nanda等[82]提出的液核生成假设,即表面先熔化随着温度逐渐向内核区域扩散。上述融化机制表明,在达到融点之前,表面预融的硼颗粒会和周围的颗粒发生烧结、团聚,阻碍氧化剂向颗粒内部扩散,从而导致点火延迟。基于机器学习势的融点研究从原子角度解释了硼颗粒在燃烧爆炸应用时受限的物理化学机制。

图10 铝、镁、硼的归一化融点[79]

2.4 碳氢燃料

碳氢燃料具有高能量密度和低成本等特性,是火箭推进剂、喷射推进剂等航空燃油的主要组成部分,其反应机制是燃烧反应动力学中重要的研究对象。Zeng等[83]将深度势模型引入燃烧领域研究,基于MN15/6-31G**级别的第一性原理数据开发机器学习势,对甲烷氧化的基本反应机理和动力学过程开展研究。分子模拟结果表明,甲醛(CH2O)和甲氧基(CH3O)是最关键的中间产物,主要反应路径为:CH4→CH3→CH3O→CH2OH→H2CO→HCO→CO。甲烷燃烧起始反应是氧气的抽氢反应,生成甲基(CH3)和HOO自由基。随着反应进行,不断生成的OH、H和HOO自由基会与甲烷发生抽氢反应,生成大量的甲基。甲基进一步被氧化为甲醛,接着失去一个氢原子形成甲酰基(HCO—),随后可以通过碰撞解离,与氧分子反应失去氢原子,从而形成一氧化碳(CO)。上述反应路径与文献结果一致,以第一性原理的精度搭建了甲烷燃烧的反应网络。

高温下碳氢燃料的不完全燃烧是碳烟前驱体多环芳香烃(Polycyclic aromatic hydrocarbons,PAH)的重要来源。碳烟的初生受燃料结构、氧气浓度、温度等因素的影响,同时存在物理和化学机制的复杂过程。Wang等[84]通过Gaussian计算得到MN15/6-31G**级别的数据集从而搭建机器学习势函数模型,并根据化学经验设计反应物,对A1～A9多环芳烃化合物(PAH)生成路径进行模拟研究。以单环和双环PAH的生成过程为例,机器学习势函数成功预测到文献中丁二烯和乙炔反应生成苯环、五元环碳碳键断裂生成萘的过程。在生成三环芳香烃化合物的模拟中,观测到的主要反应为经典的抽氢-乙炔加成(HACA)反应,即乙炔与1-萘结合生成苊并致使五元环闭合。此外,还发现一些新的反应路径,如苯基和C4H4反应生成一条长碳链,然后通过闭环和抽氢反应生成茚。这些新的反应路径对于构建详细的碳烟生成模型具有重要意义。

Chu等[85]基于DP方案,发展了PBE/DZVP精度级别的机器学习势函数,在此基础上进一步建立了纳米反应器算法以加速反应、探究密度对于反应路径的影响。纳米反应器用一个球形压头周期性地压缩和放松体系,使体系的密度发生变化,从而改变分子碰撞频率,加快反应进程。此外,高温手段使得分子动能增加,运动更为剧烈,同样可以显著加快反应进程。图11以乙炔体系多聚化反应为例,对比了常规分子动力学计算结果和纳米反应器算法加速的计算结果。

图11 3种分子动力学模拟过程中乙炔分子数量变化及物种演化[85]

其中,方案一表示纳米反应器中的反应进程,方案二和方案三分别为1 500 K和3 000 K下传统的NVT模拟。值得注意的是,纳米反应器算法可以显著加快反应进行,其反应速率是同样温度MD结果的50倍。纳米反应器中的乙炔分子数量变化与1 500 K下缩放的结果趋势相同,其反应途径相同,仅仅存在乙炔的聚合反应。这表明,增加浓度或压力只影响碰撞频率,而不影响反应路径。高温会影响体系的总能量,进而改变反应进程,因此方案三中的反应路径与前两种方案的结果并不相同,涉及到氢抽离和氢加成反应。上述结果表明纳米反应器在不影响反应路径的前提下提高了碰撞频率,为研究成核和低温反应等缓慢动力学过程的反应机理提供了一种很有前景的工具。

作者进一步采用纳米反应器方法来加速碳烟的成核过程,考虑了文献中报道的7种不同PAH分子的成核反应进程。图12展示了从纳米反应器中的不同前体收集的最终产物。对于小分子前驱体,最终生成的多环芳香烃以气态形式排列。大多数自由基会形成交联的二聚体或三聚体,并未观察到大的团簇。质量较大的前驱体分子最终会在纳米反应器的中心形成团簇。而对于大型前驱体,其最终产物呈球形堆积结构,不同层之间存在交联。上述研究表明,碳烟的最终产物形态与前驱体的大小有很大关系。随着前驱体尺寸的增加,产物颗粒从无序的团簇变成了有序的紧密堆积结构。

图12 不同前驱体在纳米反应器中生成的最终产物[85]

3 机器学习势的计算优势

原子间相互作用势的函数形式、理论假设和近似合理性决定了其精确度,进而影响分子动力学模拟结果的可靠性[86]。一般的力场假设的函数形式较为简单,如施加两体或三体截断。目前含能材料分子动力学研究常采用的ReaxFF力场,是仅对能量项进行拟合,并未考虑原子受力项,其函数形式复杂,并行效率低。机器学习势是利用其底层算法对高维函数的强大拟合能力,以一种基于数据驱动的方式近似第一性原理原子相互作用。图13总结了目前含能材料和碳氢燃料中的机器学习势模型精度。其中,Zeng等[83]开发的甲烷燃烧机器学习势是基于Gaussian计算得到MN15/6-311G**精度的第一性原理数据,其余机器学习势的数据集均为CP2K计算的PBE/DZVP级别精度。

图13 NNP势和ReaxFF反应力场精度对比

总体来说,各个体系中的机器学习势表现均优于ReaxFF力场,能量和受力的预测精度分别提高3.5倍和6.6倍(RDX、Al和碳烟体系)。最初的ReaxFF力场开发是针对碳氢燃烧并逐步拓展到CHON类含能材料研究,经过近二十年的验证和优化,其参数精度较高、力场表现较好。从图13中可以看出,经过DP-GEN迭代的RDX和碳氢的势函数模型精度都仅仅是略高于ReaxFF。以RDX为例,在能量预测方面,NNP力场相较于ReaxFF力场能量精度提高4倍,受力精度提高9倍。这是由于ReaxFF拟合过程仅考虑了能量项,并未考虑受力项,所以对于原子受力的偏差较大。值得注意的是,对于反应力场训练集和目标体系并不匹配的力场而言,ReaxFF的表现不尽人意。以硼的机器学习势为例,ReaxFF力场的能量和受力偏差分别是NNP模型的180倍和18倍。这样显著的偏差强调了力场训练集对于目标体系的适配性和移植性。此外,势能面的采样质量对于机器学习势的精度也至关重要。对于新型含能材料ICM-102,VRMD采样方案可对势能面上的高能区域进行采样,从而提升数据集的整体质量。采用VRMD采样方案使得能量预测偏差降低58%,受力预测偏差降低37%,这种辅助采样策略极大地提高了模型的整体精度。

ReaxFF力场在拟合过程中未考虑材料在极端条件下发生压缩、膨胀等力学行为。在强烈的爆炸冲击作用下,含能材料的温度会瞬间急速上升,压力升高至几十GPa,体积减少30%[87]。机器学习势可以在训练集中加入不同压缩和拉伸条件下的高温高压数据集,这对于研究含能材料的冲击压缩效应具有现实意义。图14为不同计算方法得到的RDX晶体状态方程。其中,机器学习势预测结果与第一性原理计算结果一致。

图14 RDX晶体的状态方程曲线[44]

值得注意的是,机器学习势具有一定的外推性能。在图14中,机器学习势的训练集仅考虑了0.9～1.0的压缩比,即图中的矩形部分。但是该模型能成功预测0.8～0.9区间的状态方程曲线。相比之下,反应力场由于并未考虑势能面上的压缩区域,无法复现RDX的状态方程曲线,预测的平衡结构体积偏大,在进行弛豫后将无法得到准确的晶体构型。此外,反应力场明显高估了压缩过程中的势能,在爆炸或冲击载荷的研究中可能会人为形成局部热点,加速反应的发生。因此,目前的反应力场的精度并不能满足爆炸冲击下材料性能的研究需求。Chu等[44]提出的机器学习力场更能准确描述RDX晶体在燃烧爆炸等极端环境下的理化性质和分解过程。

除了拥有更高的计算精度,机器学习势计算效率也更为优异。图15以RDX晶体为例对比了不同方法的计算效率,体系规模从100个至10万个原子。机器学习势在保持DFT精度的基础上,计算效率与体系规模呈线性关系,约为反应力场的27倍,比DFT效率提升了4个数量级。这是由于机器学习势模型计算是基于GPU硬件实现的。对于RDX晶体而言,机器学习势的计算成本与体系规模呈O(N)的趋势。而对于密度较低的气相碳反应体系,机器学习势模型的计算成本与体系规模呈O(NlogN)的趋势。机器学习势的优异表现使得用第一性精度探索具有数十万甚至数千万原子系统成为可能,可以从原子角度研究含能材料的复杂反应机制,为建立详细的含能材料反应动力学模型提供了新的机遇。

图15 密度泛函理论(DFT)、反应力场(ReaxFF)和机器学习势(NNP)的计算效率对比[44]

4 挑战与展望

分子动力学模拟能够在原子尺度上揭示含能材料的微观化学反应路径,为研究复杂体系的反应动力学机理提供了一条新的途径。分子动力学模拟的关键是力场参数,它能够描述反应的势能面结构,决定了模拟的精度和效率。从头算分子动力学准而不快,反应力场计算快而不准,机器学习势兼顾计算精度和计算效率,在含能材料燃烧与爆炸领域得到广泛应用。

机器学习算法不需要预设物理化学规律,可直接从计算数据出发,准确推断函数关系。得益于这一优点,机器学习势往往能够重现含能材料反应过程势能面变化,进而研究其热解或燃烧过程的详细反应机理。同时,数据驱动这一特性也给机器学习势的应用带来了一系列的挑战。机器学习势的开发需要大量高质量数据,且很难先验地判断数据集的完备性。虽然,数据集的规模可能很大(大于10万量子化学计算数据),但它并不足以对目标体系进行充分采样。例如,基于低温分子模拟采样搭建的数据集,其机器学习势无法准确外推到高温体系。因此,数据集需要尽可能地包括目标体系的整个势能空间。其次,训练集一般是基于第一性原理计算开发,而第一性原理计算本身也存在一定误差,受到基组大小、泛函精度等因素的限制。通过提高第一性原理计算的理论方法可以减少计算误差,然而也会极大地增加计算耗时,尤其是考虑到现有的势函数开发往往需要十万量级的第一性原理计算数据,搭建高精度数据集的成本极其昂贵。未来含能材料的机器学习势函数开发面临着如下挑战:

一是如何对极端条件下复杂反应势能面进行充分采样。此前的工作主要针对单质体系,在实际应用中含能材料的配方日益复杂,包含硝胺类材料、高分子化合物、金属材料等多种组分。随着元素种类的增加,反应空间构型的数量呈现指数级增长,需要包含更多的反应路径、中间体构型。这对如何有效地构建训练集提出了新的要求。目前存在主动学习、势能面扫描等思路。Zhang等[54]提出了DP-GEN主动学习算法,通过比较不同模型的不确定性动态生成新的数据点,从而提高模型对全局势能面的覆盖度。另一方面,研究者可以根据对化学反应的认知,设计势能面扫描算法对化学反应路径上的结构构型进行采样。Liu等[88]提出采用随机势能面行走方法(SSW),利用化学键振动对构型施加微扰,从而高效搜索复杂结构的全局势能面。Zhang等[89]提出了一种准经典轨迹法构建局部势能面的方法,可针对特定反应开发高精度势能面。增强采样方法也是一种高效采样方法,可以作为常规采样方法的重要补充。本课题组最近结合虚拟现实技术开发了一种增强采样工作流[52],在虚拟现实空间中直接对关键反应路径的势能面进行采样。和传统的分子动力学模拟直接采样相比,上述几种方法都能够提高势能面采样的效率,克服势能面上存在的高势垒,更好地描述复杂反应过程。值得注意的是,含能材料应用中常涉及极端的高温高压环境,在构建训练数据集时需要把极端条件的构型纳入考虑。

二是如何提高机器学习势训练集的理论精度。机器学习势的精度依赖于训练集选取的量子化学计算方法,特别是特定的物理化学性质或反应体系对DFT的泛函精度要求较高(如反应能垒、生成焓等),需要使用高精度的杂化甚至双杂化泛函计算方法。然而高精度的DFT计算成本随着系统规模的扩大而急剧增加,这限制了高精度DFT数据的数量,难以满足机器学习势开发的需求。Chen等[90]提出了DeePKS方法,通过机器学习方法开发DFT泛函,从而极大地缩减了高精度DFT计算的耗时,有望实现大规模的高精度DFT训练集构建。Smith等[91]提出了一种程序框架实现了小样本数据的机器学习势函数开发。他们先使用普通的DFT数据训练模型,然后使用迁移学习方法在高精度的DFT数据上对模型进行再训练,从而获得高精度的机器学习势函数。在未来开发含能材料的机器学习势函数时,结合上述方法有希望更好地构建高精度DFT数据集。