李昌霖，甘 強，馮長根，胡靖偉，朱雙飛，程年壽

（北京理工大學(xué)爆炸科學(xué)與技術(shù)國家重點實驗室，北京 100081）

1 引言

籠形多硝基化合物六硝基六氮雜異伍茲烷（簡稱HNIW 或CL-20）［1］是已應(yīng)用能量最高的單質(zhì)炸藥，由于分子中硝基空間取向和晶格堆積方式不同，常溫常壓下CL-20 主要存在α、β、γ 和ε 四種晶型，其中ε 晶型密度最大。在極端條件下CL-20 可發(fā)生復(fù)雜相變反應(yīng)，導(dǎo)致密度減小或內(nèi)部缺陷產(chǎn)生，影響其爆轟性能和安全性。實驗研究表明，β-CL-20 在高溫和高壓下發(fā)生β→ε 和β→α 相變，γ-CL-20 可由β-和γ-CL-20 轉(zhuǎn)化而成［2］，高溫下發(fā)生ε→γ［3-5］和α→γ 相變［6］，高壓下存在ε→γ［7-8］相變。由于升溫速率等實驗條件不同，已報道CL-20 相變溫度存在一定差異，且實驗方法一般只能適用于CL-20 兩相體系相變分析，難以考察晶體中分子微觀變化。2016 年張力［9］采用分子動力學(xué)方法和凝聚態(tài)優(yōu)化分子力場（COMPASS）計算表明，ε-和β-CL-20 相變溫度分別為450 K 和490 K，略高于實驗值［10］。

2001 年，van Duin 等［11］首次在分子動力學(xué)模擬方法中引入鍵級概念，提出針對碳氫化合物的反應(yīng)性力場（Reactive Force Field，ReaxFF），2011 年Liu 等［12］對該力場的長程色散作用進行低梯度校正，提出適用于含能材料的ReaxFF-lg 力場，對黑索今（RDX）、奧克托今（HMX）等含能材料計算晶胞參數(shù)誤差≤1%，已應(yīng)用于HMX 的熱膨脹和相變研究［13］，尚未應(yīng)用于CL-20高 溫 相 變 分析。Wang 等［14］采用ReaxFF-lg 力 場 研 究了ε-CL-20 熱分解機理，但常態(tài)下計算密度只有1.921 g·cm-3，與實驗值（2.035 g·cm-3）［15］偏差較大；Ren等［16］采用ReaxFF-lg力場研究β-CL-20熱分解機理，計算密 度 為1.938 g·cm-3，也與實 驗 值（1.989 g·cm-3）［17］有一定偏差。

本研究對ReaxFF-lg 力場的價角勢能截距進行了修正，計算了ε-、β-和γ-CL-20 的晶胞參數(shù)、晶格能和升華焓，采用Birch-Murnaghan 固體狀態(tài)方程擬合了ε-CL-20 的p-V 狀態(tài)曲線，驗證力場適用性。通過晶胞參數(shù)、晶體結(jié)構(gòu)變化考察了ε-、β-和γ-CL-20 高溫相變過程，分析CL-20 熱膨脹規(guī)律。

2 實驗部分

2.1 ReaxFF-lg 力場表達式及價角勢能修正

采用ReaxFF-lg 力場開展分子動力學(xué)模擬，其系統(tǒng)總能量由下列價鍵作用組成［11-12］：

式中Elp、Eover、Eunder為原子本身貢獻，Ebond、EC2為雙原子體系貢獻，Eval、Epen、Ecoa為三原子體系貢獻，Etors、Econj為四原子體系貢獻，EH-bond為氫鍵貢獻，EvdW、ECoulomb為非價鍵作用力，kJ。價角能Eangle是Eval的主要組成部分，力場將單原子與兩鍵認作鍵角體系，由兩鍵的鍵級BOa、BOb與兩鍵所成角度θ 來表達價角能的值，價角能Eangle表達式為：

式中，λ3為力場參數(shù)決定的價角參數(shù)；ka為特定角的首要力常數(shù)；BOa、BOb為兩鍵鍵級；θ 為鍵角角度，θ0價角的平衡角度，單位為（°）。當組成此角的兩鍵中的任意一鍵鍵級小于此截距時，反應(yīng)力場便會忽略此角，其缺省值為10-3。該力場中存在鍵級為10-3～10-2的弱鍵存在于同分子的兩端或不同分子之間，影響低梯度范德華力的表達。所以控制價角的鍵級截距，可消除這部分弱鍵對體系能量影響，提升該實驗方法的準確度。

2.2 模型建立與晶格參數(shù)分析

模擬計算采用大規(guī)模原子分子并行模擬程序（Large-scale Atomic/Molecular Massively Parallel Simulator，LAMMPS）［18］。從劍橋晶體數(shù)據(jù)中心提取ε-CL-20［15］、β-CL-20［19］與γ-CL-20［20］晶體數(shù)據(jù)，初始單胞內(nèi)均含有4 個分子，共144 個原子。

對初始單胞三個方向分別進行3 倍、3 倍、2 倍的周期性擴增，構(gòu)建包含72 個CL-20 分子，共2592 原子的ε-、β-和γ-CL-20 超晶胞，并將力場的價角能截距由缺省值10-3修改為10-2。采用Polak-Ribiere 共軛梯度算法［21］將各晶型CL-20 超晶胞進行能量最小化，收斂精度為4.19×10-6（kJ·mol-1）·?-1，獲得0K 時的優(yōu)化后的原子位置。隨后采用Maxwell-Boltzmann 隨機分布，給予各體系內(nèi)每個原子一個隨機初始速度。在正則（Canonical Ensemble，NVT）系綜下弛豫10 ps 后，在298 K 和0 GPa 條件下，使用等溫等壓（Isobaric-Isohermal Enermble，NPT）系綜弛豫100 ps，得到充分弛豫后結(jié)構(gòu)。NVT 與NPT 系綜所依據(jù)的運動方程均為結(jié) 合Martyna 靜 水 力 學(xué) 方 程［22］與Parrinello-Rahman應(yīng)變能方程［23］的Shinoda 運動方程［24］，并遵循Tuckerman 時間積分方案［25］。以上弛豫過程中時間步長均為0.1 fs，每一步生成一次原子近鄰列表，取最后10 ps的數(shù)據(jù)進行統(tǒng)計分析，繪制密度與體系總能量隨模擬時間的趨勢變化曲線。

2.3 升華焓

將單個分子周圍留有足夠空間可得到近似孤立分子模型，認為其與周期性邊界條件周圍分子無任何相互作用力。以ε-CL-20 為例，在NPT 系綜和700 K、0 Pa 壓力下反應(yīng)20 ps，時間步長為0.1 fs。因分子尚未發(fā)生分解，取5～20 ps 平衡后能量均值除以分子數(shù)，作為凝聚態(tài)單分子能量。創(chuàng)建與上述晶胞大小相同的空盒，僅放入一個CL-20 分子，通過周期性邊界條件預(yù)測孤立分子的能量，持續(xù)運算250 ps，取最后10 ps 軌跡的平均能量作為該分子的能量。該晶體分子晶格能（Elatt）由式（3）計算［27］：

式中，Es為單個固體分子能量，Eg為空盒子孤立分子所碎裂成原子的總能量，kJ。晶格能Elatt與升華焓成線性關(guān)系，根據(jù)Svatopluk 等［27］的含能材料晶格能和升華焓關(guān)系式Elatt=-1.2187ΔHsub+8.1804，對各晶型CL-20升華焓進行預(yù)測。

2.4 p-V 曲線擬合

為避免使用NPT 系綜時，各方向獨立控壓可能導(dǎo)致的超晶胞在高壓下的異常變形現(xiàn)象，參考低壓實驗［28］與高壓DFT 模擬結(jié)果［29］的晶胞參數(shù)值，在NVT系綜下改變盒子大小以控制壓縮率，計算298 K 時模型所受的壓力得到不同壓縮比的p-V 曲線，使用連續(xù)介 質(zhì) 中 的 三 階Birch-Murnaghan 固 態(tài) 方 程［30-31］進 行擬合：

式中，p 為壓力，GPa；B0與B′0分別為晶體的體積模量及其對壓力偏導(dǎo)數(shù)；V0/V 為晶體初始體積與在該壓力下的體積的比值。擬合三階Birch-Murnaghan 固體狀態(tài)方程，預(yù)測ε-、β-和γ-CL-20 的體積模量與其對壓力偏導(dǎo)數(shù)，與已有報道對比。

2.5 均方位移與分子質(zhì)心徑向分布函數(shù)

編寫程序計算ε-CL-20 分子的質(zhì)心，計算相對初始原子位置的均方位移與分子質(zhì)心的徑向分布函數(shù)。因衍射實驗值晶體模型與模擬優(yōu)化值晶體模型的結(jié)構(gòu)原子數(shù)量與順序均相同，因此無需考慮質(zhì)量加權(quán)。

不考慮質(zhì)量加權(quán)的均方位移（RMSD）定義如下［32］：

式中，u、v為兩原子的位置坐標，i為原子序號。使用模擬值與初始值進行比對，RMSD 越小，代表兩者結(jié)構(gòu)越加相似，用該力場表示CL-20分子結(jié)構(gòu)的置信度就越高。

同時，為驗證密度突增現(xiàn)象是否對應(yīng)晶胞內(nèi)分子分布的變化，將相變溫度對應(yīng)密度突變前后狀態(tài)的ε-CL-20 在NPT 系綜下弛豫40 ps，計算晶胞內(nèi)各個分子不同時刻的質(zhì)心坐標，對距原點分子質(zhì)心0～20 ? 范圍內(nèi)，的圓環(huán)內(nèi)分子質(zhì)心的出現(xiàn)概率進行統(tǒng)計（距離精度0.1 ?），歸一化處理后繪制分子質(zhì)心徑向分布函數(shù)（RDF）曲線。

2.6 熱膨脹及相變分析

充分弛豫得到平衡態(tài)ε-、β-和γ-CL-20 超晶胞，在0 GPa 和298～473 K 溫度下進行40 ps NPT 模擬，時間步長0.1 fs，統(tǒng)計不同溫度下晶胞參數(shù)與體積，分析各晶型CL-20 熱膨脹系數(shù)與可能存在的相變現(xiàn)象。為探究γ 晶型在高壓下的相變現(xiàn)象，在413 K、0～1 GPa下對γ-CL-20 進行計算。線膨脹系數(shù)和體積膨脹系數(shù)計算式如下：

式中，k 為線性系數(shù)；αn為膨脹系數(shù)（n=a、b、c、V），αn0分別為四個膨脹系數(shù)在常溫下對應(yīng)的初始值，K-1。

3 結(jié)果與討論

3.1 力場適用性分析

3.1.1 常態(tài)下CL-20 晶胞參數(shù)

表1 給 出 了 在298 K、0 GPa 條 件 下，ε-、β-和γ-CL-20 的晶胞參數(shù)和密度，并給出與實驗值相比［17］的RMSD 值。

表1 各晶型CL-20 晶胞的平衡參數(shù)、密度與實驗值對比Table 1 Comparison of equilibrium parameters，densities with experimental values of CL-20 cells

修正ReaxFF-lg 力場的價角勢能截距，可減小非現(xiàn)實弱鍵對分子晶體體系能量的影響，提升該力場的準確度。從表1 可見，本文得到常溫常壓下ε-、β-和γ-CL-20 晶胞參數(shù)和密度與實驗值誤差低于2%，且三種晶型CL-20 與實驗值的RMSD 值均在0.3 ? 左右，證明兩者無幾何偏差［32］，說明修正的ReaxFF-lg 力場對CL-20 的 適 用 性。Wang 等［14］采 用ReaxFF-lg 力 場 計算常溫下ε-CL-20 密度（1.921g·cm-3）與實驗值相差較大，可能是由于力場未經(jīng)價角勢能截距優(yōu)化，導(dǎo)致體系進入亞穩(wěn)定狀態(tài)。

3.1.2 晶格能與升華焓

計算得到ε-、β-和γ-CL-20 的晶格能（EL）與升華焓（Esub），如表2 所示。由表2 可見，三種晶型CL-20 晶格能與Wei 等［26］報道的晶格能基本一致，誤差在3%之內(nèi)，優(yōu)于Liu［33］的模擬結(jié)果。采用Svatopluk 關(guān)系式得到ε-、γ-和β-CL-20 的升華焓，其中ε-CL-20 的升華焓與Dorofeeva 等［34］在700 K 附 近 預(yù) 測 的 升 華 焓 相 符合。對比三種晶型CL-20 的晶格能可見，β-與γ-CL-20的晶格能大于ε-CL-20，說明ε 晶型穩(wěn)定性高于其他兩種晶型，與文獻［26］報道一致。

表2 三種晶型CL-20 晶格能與升華焓預(yù)測Table 2 Prediction of lattice energy and sublimation enthalpy of CL-20

3.1.3 ε-CL-20 的p-V 曲線

采用三階Birch-Murnaghan 固態(tài)方程［30-31］擬合了0～280 GPa 壓力范圍內(nèi)ε-CL-20 的p-V 曲線，如圖1 所示。擬合得到體積模量（B0）及其對壓力偏導(dǎo)數(shù)（B′0）見表3。

從表3 可見，在不同壓力范圍內(nèi)，擬合所得ε-CL-20 的體積模量B0逐漸增大，反映了ε-CL-20 對外界均一性壓縮的抗性增大；其對壓力偏導(dǎo)數(shù)B′0趨近于4，說明隨著壓力增大，擬合結(jié)果趨近于二階Birch-Murnaghan 固態(tài)方程。Gump 等［38］通過靜水壓測試，得到較低壓力范圍內(nèi)ε-CL-20 的體積模量B0及其對壓力偏導(dǎo)數(shù)B′0分別為（34.3±0.62）GPa 和（3.57±0.31）GPa，后續(xù)研究將數(shù)據(jù)校正為（13.6±2.0）GPa和（11.7±3.2）GPa［5］。Pu 等［35］和Sorescu 等［36-37］分別采用分子對接方法和色散力校正的密度泛函理論對ε-CL-20 的體積模量B0及其對壓力偏導(dǎo)數(shù)B′0進行計算，計算結(jié)果均與Gump 等［5］的結(jié)果較接近，本研究在0～3 GPa 范圍得到B0和B′0分別為（18.8±0.5）GPa 和（9.5±0.7）GPa，趨近于上述實驗和模擬結(jié)果。Xu［29］等模擬了0～400 GPa 范圍內(nèi)ε-CL-20 的壓縮特性，高壓數(shù)據(jù)跟本研究擬合曲線較吻合，但對三階Birch-Murnaghan 固態(tài)方程擬合度較差。

圖1 ε-CL-20 的p-V 擬合曲線Fig.1 p-V curve of ε-CL-20

如上所述，研究依據(jù)晶胞參數(shù)、RMSD、晶格能、升華焓和p-V 狀態(tài)函數(shù)等方面，說明了修正的ReaxFF-lg力場對各晶型CL-20 的適用性。上述結(jié)果說明本研究采用的修正ReaxFF-lg 反應(yīng)力場能夠可靠描述CL-20的壓縮特性。

表3 ε-CL-20 的體積模量B0及其對壓力偏導(dǎo)數(shù)B′0Table 3 Results of bulk modulus（B0）and partial derivative to pressure（B′0）of ε-CL-20

3.2 CL-20 高溫相變分析

3.2.1 晶胞參數(shù)分析

考 察 了 在298～473 K 溫 度 范 圍 內(nèi)ε-、β-和γ-CL-20 的密度變化，以及413 K 下γ-CL-20 密度隨壓力變化，如圖2 所示。對423 K 溫度下ε-CL-20 超晶胞中分子形貌變化進行了分析，如圖3 所示。

在升溫過程中，所考察三種晶型CL-20 晶胞參數(shù)與體積隨溫度上升均發(fā)生突變。從圖2a 和圖3 可見，ε-CL-20 在423 K 時發(fā)生晶胞參數(shù)突變，表現(xiàn)為密度快速下降，并略低于同溫度γ 晶型。分子排布更為疏松，CL-20 分子中部分硝基的取向從赤道轉(zhuǎn)向軸向，即轉(zhuǎn)變?yōu)棣?晶型。如圖2b 所示，將γ 晶型超晶胞在413 K下加壓，在0.5 GPa 發(fā)生晶胞密度突增，超晶胞中部分CL-20 分子的一個硝基從軸向轉(zhuǎn)向赤道，即轉(zhuǎn)變?yōu)棣啪?，?dǎo)致分子排布更加密集。Bircher 等［4］研究表明ε→γ 相變發(fā)生在347～437 K 升溫過程中；Gump 等［5］提 出398 K 時 發(fā) 生ε→γ 相變，413 K 和0.4 GPa 條件下發(fā)生γ→ε相變。本文采用價角勢能修正的ReaxFF-lg力場下得出423 K 時發(fā)生ε→γ 相變，413 K 和0.5 GPa下發(fā)生γ→ε相變，與上述兩實驗結(jié)果接近。

由圖2a 可見，β 晶型在溫度至448 K 時發(fā)生密度的突變，表現(xiàn)為分子的排布更加緊密，密度略有增大，推測發(fā)生β→ε 相變，導(dǎo)致密度增大。Russell 等［2］提出413～453 K 下發(fā)生β→ε 相變，與研究結(jié)果相符。

3.2.2 晶體結(jié)構(gòu)變化

對398 K 與423 K 兩個溫度下ε-CL-20 超晶胞進行了RDF 分析，結(jié)果見圖4。

圖2 各晶型CL-20 密度隨溫度或壓力變化曲線Fig.2 Density curve of CL-20 with temperature or pressure

圖3 423 K 時ε-CL-20 晶胞內(nèi)分子形貌變化Fig.3 Change of molecular morphology in ε-CL-20 at 423 K

圖4 398 K 與423 K 下ε-CL-20 晶體的RDF 分析Fig.4 RDF analysis of ε-CL-20 at 398 K and 423 K

觀察升溫過程中ε-CL-20 晶體內(nèi)分子變化，發(fā)現(xiàn)423 K 時ε-CL-20 晶胞內(nèi)部分分子中一個硝基從赤道位置轉(zhuǎn)向軸向，證明發(fā)生ε→γ 相變（見圖3）。RDF 曲線中第一波峰表示兩個鄰近分子質(zhì)心的距離。根據(jù)Bolotina［15］的CL-20 晶體衍射數(shù)據(jù)在相應(yīng)溫度下進行RDF 分 析，結(jié)果表明ε 與γ 晶型CL-20 晶體內(nèi)鄰近分子質(zhì)心距離分別為7.0 ? 與7.8 ?。從圖4 可見，從398 K升 溫到423 K 時，ε-CL-20 體 系 的RDF 第 一 波 峰位 置從7.3 ? 變?yōu)?.8 ?，超晶胞內(nèi)分子質(zhì)心分布發(fā)生明顯變化，也證明升溫過程中發(fā)生ε→γ 晶型轉(zhuǎn)變。

3.2.3 熱膨脹系數(shù)分析

對298～398 K 溫 度 范 圍 內(nèi)ε-、β-和γ-CL-20 晶 胞參數(shù)和體積進行線性擬合，得到各晶型CL-20 線膨脹系數(shù)與體積膨脹系數(shù)，如表4 所示。

從表4 可見，ε-CL-20 三方向熱膨脹系數(shù)差異小于2.5%，無明顯各向異性，體膨脹系數(shù)較實驗值大2.9%，計算結(jié)果與實驗結(jié)果相符，證明力場對ε-CL-20分子作用力表達較為完備。Sorescu 等［39］采用RDX 勢進行熱膨脹模擬，結(jié)果a 方向的熱膨脹系數(shù)明顯小于b、c 方向，且體膨脹系數(shù)與實驗值的誤差超過50%。

β-CL-20 籠狀結(jié)構(gòu)上的氫原子在c 方向上與相鄰分子所延展的硝基較為接近，將會在兩者之間形成氫鍵，這些氫鍵的分布是交叉的，從而抑制了該方向的熱膨脹，表現(xiàn)為c 方向的熱膨脹系數(shù)遠遠小于其余兩個方向［35］。表4 中β-CL-20 在c 方向的熱膨脹系數(shù)略微小于a、b 方向，相比于實驗值有所差異，這可能是由于力場對β-CL-20 分子間作用力的表達精度不足。而Sorescu 等［39］的模擬結(jié)果可見a 方向熱膨脹系數(shù)明顯大于b、c 方向，與實驗值不符。

表4 ε-、β-和γ-CL-20 熱膨脹系數(shù)Table 4 Coefficient of thermal expansion of ε-，β- and γ-CL-20

Pu 等［35］實驗研究表明，γ-CL-20 在b 方向的分子堆積方式按分子質(zhì)心連線可構(gòu)成六元環(huán)單元，其中垂直于b 的H—O…N 分子間作用力在同類型鍵中較弱，導(dǎo)致熱膨脹時六元環(huán)優(yōu)先沿弱鍵方向膨脹，六元環(huán)沿b 方向畸變而產(chǎn)生負膨脹現(xiàn)象。本研究的結(jié)果表明，γ-CL-20 在b 方向熱膨脹系數(shù) 顯 著 小 于a 和b 方 向，數(shù)值相差達60%，表現(xiàn)出明顯各向異性，與Pu 等的實驗結(jié)果一致，但未表現(xiàn)出負膨脹的現(xiàn)象。而Sorescu等［39］計算結(jié)果中b 方向熱膨脹系數(shù)反而有所增大，與實驗值不符。

如上所述，本研究模擬得到ε-、β-和γ-CL-20 膨脹系數(shù)與Pu 等［35］的實驗值較一致，顯著優(yōu)于Sorescu等［39］的熱膨脹計算結(jié)果。由于ReaxFF-lg 力場對氫鍵作用的表述精度有待提高，導(dǎo)致了β-和γ-CL-20 的熱膨脹各向異性差異不顯著，有待對力場進一步改進。

4 結(jié)論

對ReaxFF-lg 力場進行價角勢能截距修正，驗證了該力場對含能材料CL-20 適用性，并在該力場下對CL-20 展開了熱膨脹和相變研究，取得如下結(jié)論：

（1）計算了常溫常壓下ε-、β-和γ-CL-20 的晶胞參數(shù)、晶格能和升華焓，以及0～280 GPa 范圍內(nèi)ε-CL-20的p-V 曲線，表明修正的ReaxFF-lg 反應(yīng)力場能夠可靠描述CL-20 理化性質(zhì)。

（2）高溫相變研究表明，ε-CL-20 在423K 時發(fā)生ε →γ 相 變，β-CL-20 在448 K 時 轉(zhuǎn) 變 為ε 晶 型，而γ-CL-20 在413 K 和0.5 GPa 轉(zhuǎn)變?yōu)棣?晶型。

（3）熱膨脹系數(shù)分析表明，ε-CL-20 在升溫膨脹過程中無明顯各向異性，與實驗結(jié)果一致。而β 與γ 晶型CL-20 分別在c 方向和b 方向的熱膨脹明顯小于其他方向，表現(xiàn)出各向異性，與實驗值有一定差異，需對ReaxFF-lg 的氫鍵作用進一步修正。