几种多态含能材料构型转变的DFT研究

姚如意，张树海，苟瑞君，李连强

(1.中北大学环境与安全工程学院，山西 太原 030051；2.兵器工业卫生研究所，陕西 西安 710065)

引 言

多态含能材料晶型转变的研究对炸药的爆轰性能和安全性具有重要意义。对几种多态含能材料在熔态TNT和DNP中的溶解性研究表明[1]，RDX、HMX及CL-20在熔态TNT中的溶解度分别为3.47g /100g、0.24g /100g和5.78g /100g，在熔态DNP中的溶解度分别为12.28g /100g、2.64g /100g和7.02g /100g 。而一旦高能固相炸药在熔铸炸药载体中能溶解，则高能固相炸药在熔态的熔铸炸药载体中就可能会随着温度的变化溶解结晶，此过程还有可能造成多态含能材料晶型的转变。在对几种多态含能材料在熔态TNT和DNP中的溶解性及其结晶晶型的研究中发现[1]，CL-20在TNT、DNP中溶解后回收的晶型均由ε型变为β型。弹药中炸药发生相变时，其密度、感度、热稳定性等就会随之发生变化，导致晶体体积发生膨胀或收缩，在晶体内部形成内应力和损伤缺陷，成为潜在的热点和剪切带，影响武器弹药的使用和贮存[2]。因此多态含能材料晶型转变的研究对于熔铸炸药生产、贮存、运输和使用过程的稳定性及其相应的武器弹药安全性、可靠性研究均有重要意义[3]。

多态含能材料晶型转变的研究多为针对始态和末态稳定存在的晶型在热力学和动力学角度进行计算和解释，对多态含能材料构型转变的路径进行的计算却很少。Mrinal Ghosh等[4]采用QST2的方法计算了β-RDX→α-RDX、γ-HMX→β-HMX、ε-CL-20→β-CL-20构型转变过程。本研究采用TS方法，搜寻了β-RDX→α-RDX、γ-HMX→β-HMX、ε-CL-20→β-CL-20及β-FOX-7→α-FOX-7分子构型转变过程中的过渡态结构，确定它们的转变路径，通过研究分子构型的转变过程在一定程度上可以大致反映出晶型在转变过程中分子构型变化的一个基本情况。并通过计算吉布斯自由能随构型转变过程的变化，分析多态含能材料分子构型转变的难易顺序。

1 理论计算

1.1 过渡态搜寻

1.1.1 计算环境

多态含能材料在溶剂中之所以容易发生晶型转变，与溶剂分子和多态含能材料分子间作用力分不开。为了更有效地对多态含能材料的构型变化进行分析，本研究在溶剂丙酮中对不同分子构型间的过渡态结构进行搜寻，为了更好地表现溶剂的平均效应，计算时采用隐式溶剂环境。

1.1.2 计算方法

借助Gaussian 09[5]软件，运用密度泛函理论(DFT)，首先在密度泛函B3LYP方法和基组6-31G*下于隐式溶剂丙酮中对β-RDX与α-RDX[6-7]、γ-HMX与β-HMX[8-9]、ε-CL-20与β-CL-20[10]、α-FOX-7与β-FOX-7[11]进行结构优化；然后进行势能面柔性扫描，快速确定过渡态结构在势能面上的大致位置区间，初猜过渡态的结构；在B3LYP/6-311+G(D，P)及MO62X/6-311+G(D，P)方法和基组下，选用LTS算法，对β-RDX→α-RDX、γ-HMX→β-HMX、ε-CL-20→β-CL-20和β-FOX-7→α-FOX-7构型间的过渡态构型进行优化，直至频率分析中有且仅有一个虚频，此时初步确定是反应的过渡态结构；通过IRC计算，对过渡态结构进行验证，通过此方法搜寻出分子构型在转变过程中的所有过渡态结构。图1为优化后的α-RDX、β-RDX、β-HMX、γ-HMX、β-CL-20、ε-CL-20、β-FOX-7和α-FOX-7的分子构型图。

图1 分子构型图

1.2 吉布斯自由能计算

借助Gaussian 09软件，首先在B3LYP/6-311+G(D，P)的方法和基组下对α-RDX与β-RDX、β-HMX与γ-HMX、ε-CL-20与β-CL-20、α-FOX-7与β-FOX-7及它们各自的过渡态结构气相环境下进行优化，继续在B63LYP/6-311+G(D，P)的方法和基组下进行热力学参数计算，然后将得到的热力学校正量加到PWPB95-D3/def2-QZVP方法下计算的电子能量(E)上。B3LYP/6-311+G(D，P)方法和基组下的零点能(ZPE)、焓(ΔH)和熵(S)的校正因子[12]为0.9887、1.0102、1.0161，吉布斯自由能计算公式为G=E+ZPE+ΔH-TS。

2 结果与讨论

2.1 β-RDX→α-RDX的分子构型转变路径

2.1.1 TS计算

通过观察对比图1(a)和图1(b)，分析由β-RDX到α-RDX，主要是有一个硝基发生了偏转，在B3LYP/6-311+G(D，P)下对此键角进行势能面柔性扫描，借助Gview在β-RDX的基础上对其中一个硝基的角度不断调节作为过渡态的初猜结构，然后进行TS计算，直至出现唯一的虚频-53.05cm-1，此过渡态结构记为TS1。

2.1.2 IRC验证

对这一过渡态结构继续在相同的方法和基组B3LYP/6-311+G(D，P)下进行IRC计算，从计算结果分析TS1处于势能面的能量最高点，从最高点开始分别向曲线两侧趋向β-RDX和α-RDX结构，因此确定TS1即为β-RDX→α-RDX构型转变过程的过渡态结构。

2.1.3 转变过程

β-RDX→α-RDX的构型转变过程如图2所示。

图2 β-RDX→α-RDX的分子构型转变过程

由图2可知，β-RDX→α-RDX的构型转变过程为β-RDX→TS1→α-RDX。

2.1.4 吉布斯自由能

对β-RDX、TS1和α-RDX的电子能量(E)、零点能(ZPE)、焓(ΔH)、熵(S)和吉布斯自由能进行计算。β-RDX→α-RDX分子构型转变过程中的吉布斯自由能走势如图3所示。

图3 β-RDX→α-RDX分子构型转变过程中的吉布斯自由能走势

由图3可以看出，β-RDX→α-RDX分子构型转变时首先需要经过渡态结构TS1，这个转变需要克服的自由能能垒为5.25kJ/mol。

2.2 γ-HMX→β-HMX的分子构型转变路径

2.2.1 TS计算

通过观察对比图1(c)和图1(d)，分析γ-HMX和β-HMX之间主要是有两个硝基发生了偏转，猜测由γ-HMX→β-HMX的分子构型转变过程中可能不止一个过渡态结构。在M062X/6-31G*计算水平下，对这两个键角分别进行了势能面柔性扫描；借助Gview软件，在γ-HMX的基础上首先分别对这两个硝基的角度进行调节来作为过渡态结构的初猜；然后进行TS计算，直至频率分析中出现唯一的虚频；对这两个硝基的键角分别初猜结束后，再在现有的过渡态结构基础上改变另一个硝基的键角来作为下一个过渡态搜寻的初猜，直至TS计算后频率分析中有且仅有一个虚频，以此方法不断搜寻所有可能的过渡态结构。最终确定HMX两种分子构型间有3个过渡态结构TS1、TS2和TS3，虚频分别为-81.20、-59.52和-59.41cm-1。

2.2.2 IRC验证

对所有的过渡态结构在同样的方法和基组M062X/6-31G*下进行IRC计算，发现TS1、TS2和TS3的IRC能量曲线中能量最高点正是对应的过渡态结构TS1、TS2和TS3，且TS1的IRC能量曲线中趋于平缓的曲线两侧分别指向γ-HMX与TS2的方向，TS2的IRC能量曲线中趋于平缓的曲线两侧分别指向TS1与TS3，TS3的IRC能量曲线中趋于平缓的曲线两侧分别指向TS2与β-HMX的方向，以此验证TS1、TS2和TS3确实是γ-HMX→β-HMX分子构型转变过程中的过渡态结构。

2.2.3 转变过程

γ-HMX→β-HMX的分子构型转变过程如图4所示。

图4 γ-HMX→β-HMX的分子构型转变过程

由图4可知，γ-HMX→β-HMX的构型转变过程为γ-HMX→TS1→IN1→TS2→IN2 →TS3→β-HMX。

2.2.4 吉布斯自由能

对γ-HMX、TS1、IN1、TS2、IN2、TS3和β-HMX的电子能量(E)、零点能(ZPE)、焓(ΔH)、熵(S)和吉布斯自由能进行计算。γ-HMX→β-HMX分子构型转变过程中的吉布斯自由能走势图如图5所示。

由图5可以看出，γ-HMX→β-HMX的分子构型转变不是一步完成的，γ-HMX先形成自由能比其低的IN1结构，中间经过TS1结构，由γ-HMX到IN1需要的自由能能垒为22.21kJ/mol，接着IN1又变成比其自由能更低的IN2，由IN1越过TS2需要克服的自由能能垒为4.68kJ/mol，最后才由IN2越过TS3，克服5.93kJ/mol的自由能能垒转变成β-HMX。这一转变过程与Ostwald规则[13]刚好相符。

图5 γ-HMX→β-HMX分子构型转变过程中的吉布斯自由能走势

2.3 ε-CL-20→β-CL-20的分子构型转变路径

2.3.1 TS计算

通过观察对比图1(e)和图1(f)，初步分析ε-CL-20和β-CL-20之间主要是有两个硝基发生了偏转，猜测在ε-CL-20→β-CL-20的分子构型转变过程中不只一个过渡态结构。在密度泛函B3LYP和基组6-311+G(D，P)下，对这两个键角分别进行了势能面柔性扫描；借助Gview软件，在ε-CL-20的基础上分别先对这两个硝基的角度进行调节来作为过渡态的初猜结构；然后进行TS计算，直至出现唯一的虚频；再在现有的过渡态结构的基础上改变另一个硝基的键角来作为下一个过渡态结构搜寻的初猜，直至TS计算后频率分析中有且仅有一个虚频，以此方法不断搜寻所有的过渡态结构。最终确定ε-CL-20→β-CL-20分子构型转变过程中的过渡态结构为TS1和TS2，虚频分别为-48.55和-53.43cm-1。

2.3.2 IRC验证

在B3LYP/6-311+G(D，P)下对TS1和TS2进行IRC验证，发现TS1和TS2的IRC能量曲线中能量最高点正是对应的过渡态结构TS1和TS2，且TS1的IRC能量曲线中趋于平缓的曲线两侧分别指向ε-CL-20与TS2的方向，TS2的IRC能量曲线中趋于平缓的曲线两侧分别指向TS1与β-CL-20的方向，以此验证TS1和TS2确实是ε-CL-20→β-CL-20分子构型转变过程中的过渡态结构。

2.3.3 转变过程

ε-CL-20→β-CL-20的分子构型转变过程如图6所示。由图6可知，ε-CL-20→β-CL-20的构型转变过程为ε-CL-20→TS1→IN1→TS2→β-CL-20。

图6 ε-CL-20→β-CL-20的分子构型转变过程

2.3.4 吉布斯自由能

对ε-CL-20、TS1、IN1、TS2和β-CL-20的电子能量(E)、零点能(ZPE)、焓(ΔH)、熵(S)和吉布斯自由能进行计算。ε-CL-20→β-CL-20分子构型转变过程中的吉布斯自由能走势如图7所示。

图7 ε-CL-20→β-CL-20分子构型转变过程中的吉布斯自由能走势

由图7可知，ε-CL-20→β-CL-20的分子构型转变不是一步完成的，ε-CL-20先越过TS1，克服自由能能垒9.69kJ/mol变成IN1，最后IN1才越过TS2变成ε-CL-20，由IN1转变为TS2需要克服的自由能能垒为6.28kJ/mol。

2.4 β-FOX-7→α-FOX-7的分子构型转变路径

2.4.1 TS计算

通过观察对比图1(g)和图1(h)，初步分析β-FOX-7和α-FOX-7之间的差异主要表现在两对二面角上，在B3LYP/6-311+G(D，P)方法和基组下，对这两个键角分别进行了势能面柔性扫描；借助Gview软件，在β-FOX-7的基础上首先分别对这两个硝基的角度进行调节来作为过渡态结构的初猜；然后进行TS计算，直至频率分析中出现唯一的虚频；再在现有的过渡态结构的基础上改变另一个硝基的键角来作为下一个过渡态搜寻的初猜，接着进行TS计算。最终确定由β-FOX-7→α-FOX-7分子构型转变过程中有一个过渡态结构TS1，虚频为-101.27cm-1。

2.4.2 IRC验证

在B3LYP/6-311+G(D，P)下对TS1进行IRC验证，IRC能量曲线中能量最高点正是对应的过渡态结构TS1，且TS1的IRC能量曲线中趋于平缓的曲线两侧分别指向β-FOX-7与α-FOX-7的方向，以此验证TS1确实是β-FOX-7→α-FOX-7分子构型转变过程中的过渡态结构。

2.4.3 转变过程

β-FOX-7→α-FOX-7的分子构型转变过程如图8所示。

图8 β-FOX-7→α-FOX-7的分子构型转变过程

由图8可知，在β-FOX-7→TS1→α-FOX-7的分子构型转变过程中，原子1-2-5-6所在的二面角变化为21°→-6°→-34°，原子3-4-5-6所在的二面角变化为21°→7°→-6°。β-FOX→α-FOX的构型转变过程为β-FOX→TS1→α-FOX。

2.4.4 吉布斯自由能

对β-FOX-7、TS1和α-FOX-7的电子能量(E)、零点能(ZPE)、焓(ΔH)、熵(S)和吉布斯自由能进行计算。β-FOX-7→α-FOX-7分子构型转变过程中的吉布斯自由能走势如图9所示。

图9 β-FOX-7→α-FOX-7分子构型转变过程中的吉布斯自由能走势

由图9可以看出，由β-FOX-7到α-FOX-7需要越过的自由能能垒为10.24kJ/mol。

3 结 论

(1)β-RDX→α-RDX、γ-HMX→β-HMX、ε-CL-20→β-CL-20和β-FOX-7→α-FOX-7的分子构型转变过程分别为β-RDX→TS1→α-RDX、γ-HMX→TS1→IN1→TS2→IN2→TS3→β-HMX、ε-CL-20→TS1→IN1→TS2→β-CL-20和β-FOX-7→TS1→α-FOX-7。由这几种多态含能材料分子构型转变过程可知，分子构型转变并非一步完成的。

(2)多态含能材料分子构型间的转变，首先会由亚稳构型越过过渡态结构，然后才会进一步向亚稳结构转变。β-RDX→α-RDX、γ-HMX→β-HMX、ε-CL-20→β-CL-20和β-FOX-7→α-FOX-7需要克服的自由能能垒分别为5.25、22.21、9.69和10.24kJ/mol，转变的难度大小排序为：HMX≫ FOX-7>CL-20>RDX。