含Cl反應力場ReaxFF的評估

肖卓君,肖 斌,何 祺,王晴晴,劉 馥,劉艷潔,王 寧,李 偉,劉 軼

(1.上海大學 材料基因組工程研究院,上海 200444;2.上海大學 理學院物理系,上海 200444;3.上海電子信息職業技術學院 機械與能源工程學院,上海 201411;4.湖北航天化學技術研究所,湖北 襄陽 441003)

引 言

固體推進劑作為常見的含能材料,在軍事、民用領域以及土木工程等動力能源材料方面都有廣泛的應用[1]。高氯酸銨(AP)是復合改性雙基推進劑和復合推進劑中常用的氧化劑,被廣泛運用于固體火箭推進劑、炸藥和煙火藥的生產中[2]。因此,理解AP熱分解的微觀機理對于推進劑的制備、加工、貯存、運輸和使用等方面都具有重要的理論參考價值。大量實驗研究[3-9]為理解AP熱分解微觀機理做出了探索。隨著模擬手段的穩步發展,原子層次的計算方法已逐漸被用于對AP性質和分解過程進行深入研究。劉博等[10]用第一性原理方法計算AP晶體的能帶結構、態密度等材料物理性能, Chu等[11]使用神經網絡模型對AP分解進行模擬。

為闡明包括反應路徑、中間和最終產物在內的復雜反應規律,反應分子動力學作為一種有效研究大尺度凝聚態化學反應的模擬方法,是研究含能材料熱解和燃燒的有力手段。本研究通過比較力場和量子力學計算結果以及實驗事實,測試了兩種現有含Cl有機力場參數的有效性和適應性,并總結出了適用于模擬AP分解反應系統的力場參數優化、改進方向,為后續相關反應力場的開發奠定基礎。

1 計算方法

1.1 反應力場 (Reactive Force Field)

反應分子動力學(Reactive Molecular Dynamics, RMD)方法是一類能夠描述化學鍵光滑斷裂和形成的分子模擬方法,其核心是描述原子間相互作用的反應力場(Reactive Force Field,簡稱ReaxFF)。ReaxFF是構建在鍵級(Bond Order)勢函數基礎上的經典分子力場,如式(1)所示,其能量表達式考慮到分子內鍵長(Ebond)、鍵角(Eval)、二面角(Etors)等對能量的貢獻,以及過配位的能量矯正項(Eover)、欠配位的能量矯正項(Eunder)、鍵角能量的懲罰項(Epen)、共軛能量項(Econj)、分子間長程的范德華力(EvdWaals)和靜電相互作用項(ECoulomb)[12]。

Esystem=Ebond+Eover+Eunder+Eval+Epen+
Etors+Econj+EvdWaals+ECoulomb

(1)

基于鍵級的能量泛函設計使該力場能夠連續且平滑地描述化學鍵斷裂和形成過程中的能量變化。能量泛函參數通過擬合量子力學計算數據集獲得[13],在保證精度的同時獲得效率上的極大提高?；赗eaxFF的MD方法無需預先假定反應機制,同時考慮溫度、壓力等復雜極端外界環境因素的影響。因而RMD方法是研究含能材料熱解和燃燒反應的有力理論工具,該方法所能提供的信息是傳統量子力學/化學反應動力學方法的有益補充。

ReaxFF在含能材料領域獲得大量成功應用,包括對含能材料在機械和熱沖擊下化學反應的研究[14-16],以及2012年首次使用ReaxFF方法獲得對自燃推進劑界面發生的物理和化學過程的原子級理解[17],至2021年對自燃雙組元推進劑反應路徑和機理的分析研究[18]。截至2022年底,含C/H/O/N/Cl元素的反應力場文獻有19篇[19-37],其中有4篇文獻對Cl元素相關參數進行了優化。

出于適用程度的考量,本研究從上述4篇優化了Cl相關參數的文獻中選擇了兩篇文獻作為參考,并在文中稱之為ReaxFFC/H/O/N/Cl-2010[20]和ReaxFFC/H/O/N/Cl-2013[27]。其中ReaxFFC/H/O/N/Cl-2010力場在文獻中被用于對CuCl2水溶液中的銅配合物進行表征,ReaxFFC/H/O/N/Cl-2013力場被用于模擬Cl2和HCl氣體對鈦金屬和二氧化鈦進行蝕刻?；谝陨蟽煞N參考力場,使用Lammps軟件對典型含Cl小分子結構進行了靜態能量掃描計算,以及AP熱分解過程的反應分子動力學模擬。

1.2 量子力學計算

量子力學計算中,使用數值方法求解多體薛定諤方程組可以獲得體系的能量及其他性質,電子密度泛函理論(Density Functional Theory, DFT)通過構建粒子密度和體系基態物理性質的函數來簡化薛定諤方程組的求解問題,是一種基于量子力學的從頭算(ab-initio)理論,通常把基于密度泛函理論的計算稱為第一性原理計算(first-principles)。為合理描述AP體系,本研究選擇了適用于主族熱化學、動力學、非共價相互作用以及價態和Rydberg態電子激發能計算的M06-2X泛函[38],以及6-311++g**基組,對含Cl的典型小分子進行了第一性原理計算。共篩選了20種含Cl典型小分子(見圖1),按照元素組合類型分為Cl—H、Cl—O、Cl—C、Cl—N、Cl—Cl類以及多種元素組合類(Cl/C/H/O/N)。

圖1 20種典型含Cl小分子構型Fig.1 20 Kinds of typical configurations of small molecules containing Cl

使用Gaussian軟件對每一種典型含Cl小分子進行結構優化,并計算隨著分子中含Cl鍵長、鍵角、二面角大小改變時,分子能量相對于穩定結構時的變化,即能量掃描計算。能量掃描結果將作為力場計算結果的驗證參考,以分析兩種含Cl文獻力場的適用程度。

2 結果與討論

2.1 DFT-力場計算結果對比

針對20種典型含Cl小分子,DFT共計算了61條能量掃描曲線,總計1091個結構數據點,結果見表1。

表1 不同種類的DFT能量掃描數據集統計Table 1 Statistics of different types of DFT energy scanning datasets

對于每個DFT數據點所對應的分子結構,都分別使用ReaxFFC/H/O/N/Cl-2010和ReaxFFC/H/O/N/Cl-2013力場計算其靜態能量,并與DFT結果進行對比,結果見圖2。

圖2 鍵長—能量掃描曲線Fig.2 Bond length—energy scanning curves

對于Cl—H鍵,由圖2(c)可見,在間距1.0～2.1?內,力場ReaxFFC/H/O/N/Cl-2013的計算結果能夠很好地與DFT計算結果吻合,而力場ReaxFFC/H/O/N/Cl-2010吻合欠佳。

對于Cl—Cl、Cl—C,如圖2(a)、(e)所示,兩種力場并未出現明顯勢阱,且隨著原子間距的增加,DFT與力場計算結果偏差更加明顯。類似地,兩種文獻力場在計算Cl與O/N兩種原子鍵距對能量的影響時,也與DFT計算結果有較大偏差,見圖2(b)和(d)。以圖2(b)中ClO2分子能量隨Cl-O鍵長變化曲線為例,一個影響較大的因素在于:兩種力場均無法正確描述能量隨鍵長變化的定性關系,即力場計算結果中能量隨鍵長呈單一的負相關變化趨勢,沒有平衡位置的出現,而這顯然與常規勢能曲線不符。該誤差在后續的分子動力學模擬中則體現為:為了趨向更低的勢能值,Cl原子與O原子距離過大而無法正常成鍵,從而影響了含Cl—O鍵的相關中間體與產物的形成。這一點可以從RMD的物種分析結果獲得驗證:兩種力場模擬的反應過程中均無ClO2生成。

在研究鍵角對能量的影響時, ReaxFF反應力場也表現出與DFT方法差異較大的計算結果,大部分鍵角—能量曲線如圖3(a)所示,力場無法合理描述能量隨鍵角的變化趨勢。少數情況下如圖3(b)、(c)、(d)所示,在ReaxFFC/HO/N/Cl-2010力場仍然無法正確地描述能量隨鍵角變化的定性關系的情況下,ReaxFFC/HO/N/Cl-2013力場計算結果與DFT結果更加接近。

圖3 鍵角—能量掃描曲線Fig.3 Bond angle—energy scanning curves

二面角-能量掃描曲線(見圖4)與鍵角-能量變化曲線類似,大多數結果中,兩種力場均無法定性描述曲線趨勢,如圖4(a)。但值得一提的是,如圖4(b)、(c)所示,對于分子C2H5ClO和CH3Cl中的二面角C—Cl—C—H,兩種力場的計算結果都與DFT值吻合較好。

圖4 二面角—能量掃描曲線Fig.4 Torsion angle—energy scanning curves

結合圖5(d)和(e)小分子的ADCH電荷計算結果,從以上對比中還可發現,與DFT值吻合較好的力場計算結果,均是含低價態Cl小分子的掃描曲線。其中包括HCl分子的H—Cl鍵長掃描、C2H5ClO和CH3Cl分子的C—Cl—C—H二面角掃描。相反地,以圖5(a)、(b)、(c)中AP分子的各類能量掃描為例,含高價Cl的小分子,包括AP、ClO2等分子的各項掃描結果,與DFT計算結果對比均出現了不同程度的誤差。

圖5 AP小分子中:(a)鍵長-能量掃描曲線;(b)二面角-能量掃描曲線;(c)鍵角-能量掃描曲線;(d)HCl/C2H5ClO/NH4ClO4/CH3Cl中各原子的ADCH電荷值;(e)HCl/C2H5ClO/NH4ClO4/CH3Cl中各原子ADCH電荷與ReaxFF力場計算電荷值Fig.5 In AP small molecules, (a) Bond length -Energy scanning curves; (b) Torsion angle-energy scanning curves; (c) Bond angle-energy scanning curve; (d) ADCH charge values of each atom in HCl/C2H5ClO/NH4ClO4/CH3Cl; (e) ADCH charge values and the results from ReaxFF of each atom in HCl/C2H5ClO/NH4ClO4/CH3Cl

2.2 反應分子動力學模擬

實驗結果分析認為,AP在常溫常壓下以晶體形式存在,但在分解初期存在離解與升華現象[3],本研究的分子動力學模擬旨在探索力場的適用程度,分析力場是否能對相關產物的形成進行合理描述,故采用無序隨機的小分子排布方式進行建模,如圖6所示。

圖6 (a)RMD中AP小分子無序排布模型;(b)AP分子Fig.6 (a) Model of disordered arrangement of AP small molecules in RMD (b) AP molecule

模擬中總計400個AP小分子,取實驗密度1.95g/mL作為初始密度,使用ReaxFFC/H/O/N/Cl-2010和ReaxFFC/H/O/N/Cl-2013力場參數,積分步長0.1fs,NVT系綜下分別加熱300ps至1500,2000和2500K,并分別于目標溫度下NVT保溫500ps,分析不同溫度下兩種反應力場模擬AP分解時生成的產物。

早期對AP分解產物的研究發現,熱分解的產物成分取決于溫度[3, 7, 9, 39-40],分為低溫和高溫兩個階段,兩個階段的分解機理不同?；诒狙芯磕MAP分解的升溫范圍,認為反應屬于高溫階段。在溫度高于380℃時,AP分解反應方程式可參考:

此外,高溫階段還檢測到了NH3、ClO2、HCl、NOCl、N2等多種產物。

O2作為AP高溫分解階段的主要生成物之一[8],在兩種力場模擬的產物中均有體現,如圖7(a)和(b)所示。同樣地,兩種力場模擬結果均顯示了H2O的生成,如圖7(c)和(d)所示,與實驗結果一致。

圖7 ReaxFFC/H/O/N/Cl-2010力場MD模擬中,(a)O2、(c)H2O分子數量隨時間的變化曲線;ReaxFFC/H/O/N/Cl-2013力場MD模擬中,(b)O2、(d)H2O分子數量隨時間的變化曲線 Fig.7 In MD with ReaxFFC/H/O/N/Cl-2010, curves of(a) O2, (c) H2O as a function of time. In MD with ReaxFFC/H/O/N/Cl-2013, curves of (b) O2, (d) H2O as a function of time

根據文獻[41]由實驗結果推測的AP高溫反應過程,溫度對產物中O2的含量并無明顯影響,ReaxFFC/H/O/N/Cl-2013力場的模擬結果更符合這一推論。

圖8 ReaxFFC/H/O/N/Cl-2010力場和ReaxFFC/H/O/N/Cl-2013力場MD模擬中NH3和N2分子數量隨時間的變化曲線Fig.8 In MD with ReaxFFC/H/O/N/Cl-2010 and ReaxFFC/H/O/N/Cl-2013, curves of NH3 and N2 as a function of time

圖9為N2O和NO2中鍵長—能量掃描曲線的DFT和力場計算結果對比。在力場模擬產物中,大量N原子以N2的形式出現,如圖8(c)、(d)所示,而AP高溫下快速熱分解實驗中[9],N原子最終以NO2、NOCl、N2O以及少量NO的形式出現在產物中,其中NO2、N2O、NOCl這三類分子均未在力場模擬結果中體現。從圖9可以推測,與Cl—O鍵類似,力場對N—O成鍵行為的描述存在不足,導致N、O原子在分子動力學模擬過程中難以成鍵。

圖9 N2O、NO2中N—O鍵長—能量掃描曲線Fig.9 Bond-Energy scanning curve of N—O in N2O, NO2

圖10分別為兩種力場模擬中,HCl和Cl數量隨時間的變化曲線。

圖10 ReaxFFC/H/O/N/Cl-2010力場和ReaxFFC/H/O/N/Cl-2013力場模擬中HCl和Cl數量隨時間的變化曲線Fig.10 In MD with ReaxFFC/H/O/N/Cl-2010 and ReaxFFC/H/O/N/Cl-2013, curves of HCl and Cl as a function of time

對于含Cl分解產物,兩種力場均模擬出HCl分子的生成,且HCl的含量隨著溫度升高而增大,見圖10(a)和(b),這一點與高溫階段實驗事實相符。然而,實驗顯示[43],高溫情況下,Cl原子主要以Cl2的形式釋放,HCl分子作為產物并非主要類型。但在RMD模擬中,Cl原子主要以HCl的形式出現在產物里。這一模擬結果很大程度上源于力場文件中H—Cl參數能夠較為準確地描述H原子和Cl原子的成鍵情況,而Cl—Cl描述欠佳,從圖2(a)Cl—Cl鍵長的能量掃描結果對比圖也能夠看出,兩種參考力場下Cl原子均傾向于遠離而非成鍵。

除HCl分子,其他大多數Cl原子還以游離的單個原子狀態而非Cl2單質形式出現在力場模擬產物中,見圖10(c)、(d),這與力場中Cl—Cl鍵長參數的局限性有很大關聯。其次,在AP分解過程中,典型的含高價Cl中間體并未在力場計算結果中觀察到,包括ClO2、ClO3、HClO4等。說明現有的Cl元素相關參數對原子間作用力的描述仍不夠準確,還存在較大可供優化的空間。

AP具有明顯區別于其他推進劑含能材料的性質。首先,AP分子中Cl元素的化合價為+7價,具有極強的氧化性。大多數含能材料僅由C/H/O/N元素組成,高價態Cl的存在使得AP與其他含能材料在成分上有顯著差異。其次,AP分子包含H/O/N/Cl四種元素,如果考慮這些元素所有可能存在的氧化態,則AP的分解可能存在很多的過程。多種因素影響使得對AP分解過程進行分子動力學模擬時,從其他研究中遷移而來的力場可能存在無法合理描述在快速熱解反應過程中高價態Cl元素與其他有機元素的耦合反應行為問題。同時,專門針對含能材料的可用含Cl力場目前仍未被開發出來[44]。為指導開發含Cl有機反應力場,本研究使用兩種參考力場,計算了典型含Cl小分子的能量掃描曲線,以及AP分解的RMD模擬,并將力場計算結果分別與DFT結果、實驗結果進行對比評估。

分析比較能量掃描曲線和DFT量子力學計算結果,發現 ReaxFFC/H/O/N/Cl-2013力場的H—Cl參數表現更好,對含H、Cl的小分子(如HCl、CH3Cl)的成鍵斷鍵情況,以及部分二面角對能量的影響有較為合理的預測能力。然而,對于含高價Cl小分子,兩種力場均無法恰當地描述原子間距以及價角、二面角對能量的貢獻。其次,從分子動力學模擬分解產物和實驗觀測結果的對比分析得出,現有的含Cl有機力場對AP分解反應系統中非含Cl小分子(例如H2O、N2、NH3等)的行為有一定的復現能力,其中ReaxFFC/H/O/N/Cl-2013力場對HCl、部分非含Cl分解產物的模擬表現更加貼合實驗事實,反之,對于含高價Cl的小分子產物,以及含N/O元素的產物,由于力場對Cl—C/O/N原子之間、N—O原子之間缺乏適度的成鍵描述,使得分解反應模擬中難以形成含高價Cl小分子、NO2、N2O等產物,導致對實驗結果重要部分的復刻缺失。

此外,在兩種參考力場的分子動力學模擬過程中均觀察到,AP小分子群在加熱過程中就開始了反應,推測原因在于兩種力場的參數均未添加色散修正項(即—lg修正項[45]),導致加熱過程中NVT壓強過高,促使反應的發生。

3 結 論

(1)基于量子力學密度泛函理論計算,建立了含Cl有機小分子的數據集,并對現有文獻中的反應力場的表現進行了評估研究。

(2)本研究的評估結果對于使用ReaxFF反應力場模擬AP分解提供了重要警示:目前文獻中的ReaxFF力場雖然含有C/H/O/N/Cl參數,但對于較為特殊的AP體系,本工作證明了Cl參數的遷移性并不理想,與密度泛函計算結果相比誤差較大。主要誤差來源于能量表達式中分子內鍵長、鍵角和二面角等因素。因此,優化現有力場參數,特別是鍵長、鍵角和二面角的能量項相關參數,發展適用于高價Cl元素和氮氧化物NOx體系的力場參數,成為使用ReaxFF反應力場模擬AP復雜分解反應的迫切需求。本工作為未來開發更準確的含Cl反應力場提供了前期基礎和建議方向。