幾種多態(tài)含能材料構(gòu)型轉(zhuǎn)變的DFT研究

姚如意，張樹海，茍瑞君，李連強

(1.中北大學(xué)環(huán)境與安全工程學(xué)院，山西 太原 030051；2.兵器工業(yè)衛(wèi)生研究所，陜西 西安 710065)

引 言

多態(tài)含能材料晶型轉(zhuǎn)變的研究對炸藥的爆轟性能和安全性具有重要意義。對幾種多態(tài)含能材料在熔態(tài)TNT和DNP中的溶解性研究表明[1]，RDX、HMX及CL-20在熔態(tài)TNT中的溶解度分別為3.47g /100g、0.24g /100g和5.78g /100g，在熔態(tài)DNP中的溶解度分別為12.28g /100g、2.64g /100g和7.02g /100g 。而一旦高能固相炸藥在熔鑄炸藥載體中能溶解，則高能固相炸藥在熔態(tài)的熔鑄炸藥載體中就可能會隨著溫度的變化溶解結(jié)晶，此過程還有可能造成多態(tài)含能材料晶型的轉(zhuǎn)變。在對幾種多態(tài)含能材料在熔態(tài)TNT和DNP中的溶解性及其結(jié)晶晶型的研究中發(fā)現(xiàn)[1]，CL-20在TNT、DNP中溶解后回收的晶型均由ε型變?yōu)棣滦?。彈藥中炸藥發(fā)生相變時，其密度、感度、熱穩(wěn)定性等就會隨之發(fā)生變化，導(dǎo)致晶體體積發(fā)生膨脹或收縮，在晶體內(nèi)部形成內(nèi)應(yīng)力和損傷缺陷，成為潛在的熱點和剪切帶，影響武器彈藥的使用和貯存[2]。因此多態(tài)含能材料晶型轉(zhuǎn)變的研究對于熔鑄炸藥生產(chǎn)、貯存、運輸和使用過程的穩(wěn)定性及其相應(yīng)的武器彈藥安全性、可靠性研究均有重要意義[3]。

多態(tài)含能材料晶型轉(zhuǎn)變的研究多為針對始態(tài)和末態(tài)穩(wěn)定存在的晶型在熱力學(xué)和動力學(xué)角度進行計算和解釋，對多態(tài)含能材料構(gòu)型轉(zhuǎn)變的路徑進行的計算卻很少。Mrinal Ghosh等[4]采用QST2的方法計算了β-RDX→α-RDX、γ-HMX→β-HMX、ε-CL-20→β-CL-20構(gòu)型轉(zhuǎn)變過程。本研究采用TS方法，搜尋了β-RDX→α-RDX、γ-HMX→β-HMX、ε-CL-20→β-CL-20及β-FOX-7→α-FOX-7分子構(gòu)型轉(zhuǎn)變過程中的過渡態(tài)結(jié)構(gòu)，確定它們的轉(zhuǎn)變路徑，通過研究分子構(gòu)型的轉(zhuǎn)變過程在一定程度上可以大致反映出晶型在轉(zhuǎn)變過程中分子構(gòu)型變化的一個基本情況。并通過計算吉布斯自由能隨構(gòu)型轉(zhuǎn)變過程的變化，分析多態(tài)含能材料分子構(gòu)型轉(zhuǎn)變的難易順序。

1 理論計算

1.1 過渡態(tài)搜尋

1.1.1 計算環(huán)境

多態(tài)含能材料在溶劑中之所以容易發(fā)生晶型轉(zhuǎn)變，與溶劑分子和多態(tài)含能材料分子間作用力分不開。為了更有效地對多態(tài)含能材料的構(gòu)型變化進行分析，本研究在溶劑丙酮中對不同分子構(gòu)型間的過渡態(tài)結(jié)構(gòu)進行搜尋，為了更好地表現(xiàn)溶劑的平均效應(yīng)，計算時采用隱式溶劑環(huán)境。

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]進行結(jié)構(gòu)優(yōu)化；然后進行勢能面柔性掃描，快速確定過渡態(tài)結(jié)構(gòu)在勢能面上的大致位置區(qū)間，初猜過渡態(tài)的結(jié)構(gòu)；在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構(gòu)型間的過渡態(tài)構(gòu)型進行優(yōu)化，直至頻率分析中有且僅有一個虛頻，此時初步確定是反應(yīng)的過渡態(tài)結(jié)構(gòu)；通過IRC計算，對過渡態(tài)結(jié)構(gòu)進行驗證，通過此方法搜尋出分子構(gòu)型在轉(zhuǎn)變過程中的所有過渡態(tài)結(jié)構(gòu)。圖1為優(yōu)化后的α-RDX、β-RDX、β-HMX、γ-HMX、β-CL-20、ε-CL-20、β-FOX-7和α-FOX-7的分子構(gòu)型圖。

圖1 分子構(gòu)型圖

1.2 吉布斯自由能計算

借助Gaussian 09軟件，首先在B3LYP/6-311+G(D，P)的方法和基組下對α-RDX與β-RDX、β-HMX與γ-HMX、ε-CL-20與β-CL-20、α-FOX-7與β-FOX-7及它們各自的過渡態(tài)結(jié)構(gòu)氣相環(huán)境下進行優(yōu)化，繼續(xù)在B63LYP/6-311+G(D，P)的方法和基組下進行熱力學(xué)參數(shù)計算，然后將得到的熱力學(xué)校正量加到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 結(jié)果與討論

2.1 β-RDX→α-RDX的分子構(gòu)型轉(zhuǎn)變路徑

2.1.1 TS計算

通過觀察對比圖1(a)和圖1(b)，分析由β-RDX到α-RDX，主要是有一個硝基發(fā)生了偏轉(zhuǎn)，在B3LYP/6-311+G(D，P)下對此鍵角進行勢能面柔性掃描，借助Gview在β-RDX的基礎(chǔ)上對其中一個硝基的角度不斷調(diào)節(jié)作為過渡態(tài)的初猜結(jié)構(gòu)，然后進行TS計算，直至出現(xiàn)唯一的虛頻-53.05cm-1，此過渡態(tài)結(jié)構(gòu)記為TS1。

2.1.2 IRC驗證

對這一過渡態(tài)結(jié)構(gòu)繼續(xù)在相同的方法和基組B3LYP/6-311+G(D，P)下進行IRC計算，從計算結(jié)果分析TS1處于勢能面的能量最高點，從最高點開始分別向曲線兩側(cè)趨向β-RDX和α-RDX結(jié)構(gòu)，因此確定TS1即為β-RDX→α-RDX構(gòu)型轉(zhuǎn)變過程的過渡態(tài)結(jié)構(gòu)。

2.1.3 轉(zhuǎn)變過程

β-RDX→α-RDX的構(gòu)型轉(zhuǎn)變過程如圖2所示。

圖2 β-RDX→α-RDX的分子構(gòu)型轉(zhuǎn)變過程

由圖2可知，β-RDX→α-RDX的構(gòu)型轉(zhuǎn)變過程為β-RDX→TS1→α-RDX。

2.1.4 吉布斯自由能

對β-RDX、TS1和α-RDX的電子能量(E)、零點能(ZPE)、焓(ΔH)、熵(S)和吉布斯自由能進行計算。β-RDX→α-RDX分子構(gòu)型轉(zhuǎn)變過程中的吉布斯自由能走勢如圖3所示。

圖3 β-RDX→α-RDX分子構(gòu)型轉(zhuǎn)變過程中的吉布斯自由能走勢

由圖3可以看出，β-RDX→α-RDX分子構(gòu)型轉(zhuǎn)變時首先需要經(jīng)過渡態(tài)結(jié)構(gòu)TS1，這個轉(zhuǎn)變需要克服的自由能能壘為5.25kJ/mol。

2.2 γ-HMX→β-HMX的分子構(gòu)型轉(zhuǎn)變路徑

2.2.1 TS計算

通過觀察對比圖1(c)和圖1(d)，分析γ-HMX和β-HMX之間主要是有兩個硝基發(fā)生了偏轉(zhuǎn)，猜測由γ-HMX→β-HMX的分子構(gòu)型轉(zhuǎn)變過程中可能不止一個過渡態(tài)結(jié)構(gòu)。在M062X/6-31G*計算水平下，對這兩個鍵角分別進行了勢能面柔性掃描；借助Gview軟件，在γ-HMX的基礎(chǔ)上首先分別對這兩個硝基的角度進行調(diào)節(jié)來作為過渡態(tài)結(jié)構(gòu)的初猜；然后進行TS計算，直至頻率分析中出現(xiàn)唯一的虛頻；對這兩個硝基的鍵角分別初猜結(jié)束后，再在現(xiàn)有的過渡態(tài)結(jié)構(gòu)基礎(chǔ)上改變另一個硝基的鍵角來作為下一個過渡態(tài)搜尋的初猜，直至TS計算后頻率分析中有且僅有一個虛頻，以此方法不斷搜尋所有可能的過渡態(tài)結(jié)構(gòu)。最終確定HMX兩種分子構(gòu)型間有3個過渡態(tài)結(jié)構(gòu)TS1、TS2和TS3，虛頻分別為-81.20、-59.52和-59.41cm-1。

2.2.2 IRC驗證

對所有的過渡態(tài)結(jié)構(gòu)在同樣的方法和基組M062X/6-31G*下進行IRC計算，發(fā)現(xiàn)TS1、TS2和TS3的IRC能量曲線中能量最高點正是對應(yīng)的過渡態(tài)結(jié)構(gòu)TS1、TS2和TS3，且TS1的IRC能量曲線中趨于平緩的曲線兩側(cè)分別指向γ-HMX與TS2的方向，TS2的IRC能量曲線中趨于平緩的曲線兩側(cè)分別指向TS1與TS3，TS3的IRC能量曲線中趨于平緩的曲線兩側(cè)分別指向TS2與β-HMX的方向，以此驗證TS1、TS2和TS3確實是γ-HMX→β-HMX分子構(gòu)型轉(zhuǎn)變過程中的過渡態(tài)結(jié)構(gòu)。

2.2.3 轉(zhuǎn)變過程

γ-HMX→β-HMX的分子構(gòu)型轉(zhuǎn)變過程如圖4所示。

圖4 γ-HMX→β-HMX的分子構(gòu)型轉(zhuǎn)變過程

由圖4可知，γ-HMX→β-HMX的構(gòu)型轉(zhuǎn)變過程為γ-HMX→TS1→IN1→TS2→IN2 →TS3→β-HMX。

2.2.4 吉布斯自由能

對γ-HMX、TS1、IN1、TS2、IN2、TS3和β-HMX的電子能量(E)、零點能(ZPE)、焓(ΔH)、熵(S)和吉布斯自由能進行計算。γ-HMX→β-HMX分子構(gòu)型轉(zhuǎn)變過程中的吉布斯自由能走勢圖如圖5所示。

由圖5可以看出，γ-HMX→β-HMX的分子構(gòu)型轉(zhuǎn)變不是一步完成的，γ-HMX先形成自由能比其低的IN1結(jié)構(gòu)，中間經(jīng)過TS1結(jié)構(gòu)，由γ-HMX到IN1需要的自由能能壘為22.21kJ/mol，接著IN1又變成比其自由能更低的IN2，由IN1越過TS2需要克服的自由能能壘為4.68kJ/mol，最后才由IN2越過TS3，克服5.93kJ/mol的自由能能壘轉(zhuǎn)變成β-HMX。這一轉(zhuǎn)變過程與Ostwald規(guī)則[13]剛好相符。

圖5 γ-HMX→β-HMX分子構(gòu)型轉(zhuǎn)變過程中的吉布斯自由能走勢

2.3 ε-CL-20→β-CL-20的分子構(gòu)型轉(zhuǎn)變路徑

2.3.1 TS計算

通過觀察對比圖1(e)和圖1(f)，初步分析ε-CL-20和β-CL-20之間主要是有兩個硝基發(fā)生了偏轉(zhuǎn)，猜測在ε-CL-20→β-CL-20的分子構(gòu)型轉(zhuǎn)變過程中不只一個過渡態(tài)結(jié)構(gòu)。在密度泛函B3LYP和基組6-311+G(D，P)下，對這兩個鍵角分別進行了勢能面柔性掃描；借助Gview軟件，在ε-CL-20的基礎(chǔ)上分別先對這兩個硝基的角度進行調(diào)節(jié)來作為過渡態(tài)的初猜結(jié)構(gòu)；然后進行TS計算，直至出現(xiàn)唯一的虛頻；再在現(xiàn)有的過渡態(tài)結(jié)構(gòu)的基礎(chǔ)上改變另一個硝基的鍵角來作為下一個過渡態(tài)結(jié)構(gòu)搜尋的初猜，直至TS計算后頻率分析中有且僅有一個虛頻，以此方法不斷搜尋所有的過渡態(tài)結(jié)構(gòu)。最終確定ε-CL-20→β-CL-20分子構(gòu)型轉(zhuǎn)變過程中的過渡態(tài)結(jié)構(gòu)為TS1和TS2，虛頻分別為-48.55和-53.43cm-1。

2.3.2 IRC驗證

在B3LYP/6-311+G(D，P)下對TS1和TS2進行IRC驗證，發(fā)現(xiàn)TS1和TS2的IRC能量曲線中能量最高點正是對應(yīng)的過渡態(tài)結(jié)構(gòu)TS1和TS2，且TS1的IRC能量曲線中趨于平緩的曲線兩側(cè)分別指向ε-CL-20與TS2的方向，TS2的IRC能量曲線中趨于平緩的曲線兩側(cè)分別指向TS1與β-CL-20的方向，以此驗證TS1和TS2確實是ε-CL-20→β-CL-20分子構(gòu)型轉(zhuǎn)變過程中的過渡態(tài)結(jié)構(gòu)。

2.3.3 轉(zhuǎn)變過程

ε-CL-20→β-CL-20的分子構(gòu)型轉(zhuǎn)變過程如圖6所示。由圖6可知，ε-CL-20→β-CL-20的構(gòu)型轉(zhuǎn)變過程為ε-CL-20→TS1→IN1→TS2→β-CL-20。

圖6 ε-CL-20→β-CL-20的分子構(gòu)型轉(zhuǎn)變過程

2.3.4 吉布斯自由能

對ε-CL-20、TS1、IN1、TS2和β-CL-20的電子能量(E)、零點能(ZPE)、焓(ΔH)、熵(S)和吉布斯自由能進行計算。ε-CL-20→β-CL-20分子構(gòu)型轉(zhuǎn)變過程中的吉布斯自由能走勢如圖7所示。

圖7 ε-CL-20→β-CL-20分子構(gòu)型轉(zhuǎn)變過程中的吉布斯自由能走勢

由圖7可知，ε-CL-20→β-CL-20的分子構(gòu)型轉(zhuǎn)變不是一步完成的，ε-CL-20先越過TS1，克服自由能能壘9.69kJ/mol變成IN1，最后IN1才越過TS2變成ε-CL-20，由IN1轉(zhuǎn)變?yōu)門S2需要克服的自由能能壘為6.28kJ/mol。

2.4 β-FOX-7→α-FOX-7的分子構(gòu)型轉(zhuǎn)變路徑

2.4.1 TS計算

通過觀察對比圖1(g)和圖1(h)，初步分析β-FOX-7和α-FOX-7之間的差異主要表現(xiàn)在兩對二面角上，在B3LYP/6-311+G(D，P)方法和基組下，對這兩個鍵角分別進行了勢能面柔性掃描；借助Gview軟件，在β-FOX-7的基礎(chǔ)上首先分別對這兩個硝基的角度進行調(diào)節(jié)來作為過渡態(tài)結(jié)構(gòu)的初猜；然后進行TS計算，直至頻率分析中出現(xiàn)唯一的虛頻；再在現(xiàn)有的過渡態(tài)結(jié)構(gòu)的基礎(chǔ)上改變另一個硝基的鍵角來作為下一個過渡態(tài)搜尋的初猜，接著進行TS計算。最終確定由β-FOX-7→α-FOX-7分子構(gòu)型轉(zhuǎn)變過程中有一個過渡態(tài)結(jié)構(gòu)TS1，虛頻為-101.27cm-1。

2.4.2 IRC驗證

在B3LYP/6-311+G(D，P)下對TS1進行IRC驗證，IRC能量曲線中能量最高點正是對應(yīng)的過渡態(tài)結(jié)構(gòu)TS1，且TS1的IRC能量曲線中趨于平緩的曲線兩側(cè)分別指向β-FOX-7與α-FOX-7的方向，以此驗證TS1確實是β-FOX-7→α-FOX-7分子構(gòu)型轉(zhuǎn)變過程中的過渡態(tài)結(jié)構(gòu)。

2.4.3 轉(zhuǎn)變過程

β-FOX-7→α-FOX-7的分子構(gòu)型轉(zhuǎn)變過程如圖8所示。

圖8 β-FOX-7→α-FOX-7的分子構(gòu)型轉(zhuǎn)變過程

由圖8可知，在β-FOX-7→TS1→α-FOX-7的分子構(gòu)型轉(zhuǎn)變過程中，原子1-2-5-6所在的二面角變化為21°→-6°→-34°，原子3-4-5-6所在的二面角變化為21°→7°→-6°。β-FOX→α-FOX的構(gòu)型轉(zhuǎn)變過程為β-FOX→TS1→α-FOX。

2.4.4 吉布斯自由能

對β-FOX-7、TS1和α-FOX-7的電子能量(E)、零點能(ZPE)、焓(ΔH)、熵(S)和吉布斯自由能進行計算。β-FOX-7→α-FOX-7分子構(gòu)型轉(zhuǎn)變過程中的吉布斯自由能走勢如圖9所示。

圖9 β-FOX-7→α-FOX-7分子構(gòu)型轉(zhuǎn)變過程中的吉布斯自由能走勢

由圖9可以看出，由β-FOX-7到α-FOX-7需要越過的自由能能壘為10.24kJ/mol。

3 結(jié) 論

(1)β-RDX→α-RDX、γ-HMX→β-HMX、ε-CL-20→β-CL-20和β-FOX-7→α-FOX-7的分子構(gòu)型轉(zhuǎn)變過程分別為β-RDX→TS1→α-RDX、γ-HMX→TS1→IN1→TS2→IN2→TS3→β-HMX、ε-CL-20→TS1→IN1→TS2→β-CL-20和β-FOX-7→TS1→α-FOX-7。由這幾種多態(tài)含能材料分子構(gòu)型轉(zhuǎn)變過程可知，分子構(gòu)型轉(zhuǎn)變并非一步完成的。

(2)多態(tài)含能材料分子構(gòu)型間的轉(zhuǎn)變，首先會由亞穩(wěn)構(gòu)型越過過渡態(tài)結(jié)構(gòu)，然后才會進一步向亞穩(wěn)結(jié)構(gòu)轉(zhuǎn)變。β-RDX→α-RDX、γ-HMX→β-HMX、ε-CL-20→β-CL-20和β-FOX-7→α-FOX-7需要克服的自由能能壘分別為5.25、22.21、9.69和10.24kJ/mol，轉(zhuǎn)變的難度大小排序為：HMX? FOX-7>CL-20>RDX。