ReaxFF MD局部區域反應追蹤與物理性質可視化分析

唐鈺杰，鄭默，任春醒，李曉霞,＊，郭力

1中國科學院過程工程研究所多相復雜系統國家重點實驗室，北京 100190

2中國科學院大學化學工程學院，北京 100049

3中國科學院綠色制造創新研究院，北京 100190

關鍵字：ReaxFF MD；區域反應追蹤；物理性質；可視化

1 引言

反應分子動力學是將反應力場 ReaxFF(Reactive Force Field)和分子動力學MD (Molecular Dynamics)結合的分子模擬方法。基于Adri van Duin等1提出的反應力場ReaxFF，ReaxFF MD可較好地模擬復雜分子反應體系中鍵的生成和斷裂隨時間的變化，從而模擬化學反應隨時間的演化。目前，國際上主流的分子模擬方法包括分子動力學和量子力學方法。經典分子動力學方法基于牛頓力學主要描述體系的物理性質，可模擬大規模分子體系(～100000 – ～1000000原子)隨時間的動態演化；但由于原子的連接關系和電荷保持不變，無法描述化學反應。量子力學基于薛定諤方程，可精確描述體系中電子運動狀態，從而描述可能的反應路徑；但其計算代價高昂，能夠應用體系的原子規模相對較小(～100原子)。為此，基于量子力學的分子動力學模擬可考察體系鍵生成和斷裂的動態過程，但因計算能力的限制，難以應用于復雜分子體系。ReaxFF MD方法基于鍵級描述力勢中所有與成斷鍵相關的能量項，包括鍵長、鍵角、扭轉二面角、過配位和配位不足校正、氫鍵作用、其他相互作用的校正如孤對電子能、三體共軛、四體共軛等；采用基于Taper校正的修正Morse勢描述范德華非鍵作用；采用原子點電荷描述庫侖靜電作用，并且利用電負性平衡算法動態更新每個分子動力學時間步的原子電荷。由于無需預設反應路徑，可平滑描述多分子體系化學反應隨時間的演化過程，因此，ReaxFF MD方法是一種有潛力模擬和揭示復雜分子體系化學反應機理的新方法。

分析ReaxFF MD模擬結果中所蘊含的化學反應信息，對利用ReaxFF MD模擬方法認識復雜體系的反應機理極為關鍵。目前，國際上集成了ReaxFF MD模擬計算程序的主流平臺較多，有LAMMPS2、AMS (原名是ADF)3、Materials Studio(MS)，但專門分析ReaxFF MD模擬結果中化學反應的程序系統仍然十分缺乏。例如，MS中GULP模塊的ReaxFF 6.04雖然可以進行反應體系的模擬，但并未提供專門針對ReaxFF MD化學反應信息的分析工具，其自帶的工具主要用于分析經典分子動力學模擬的物理過程，只能觀察ReaxFF MD模擬得到的分子體系軌跡變化，不能獲得化學反應信息。從發表的ReaxFF MD模擬應用文章可知，ReaxFF MD模擬結果化學反應工具主要為LAMMPS平臺開源的reax/c模塊5，其反應信息分析結果為基于分子式的分子數量隨時間的演化，既不能直接給出反應細節(特別是反應位點的信息)，也不能區分同分異構體。

隨著ReaxFF力場和模擬應用的快速發展，在商業化的分子動力學模擬平臺中集成和發展專門用于分析ReaxFF MD模擬結果的程序在最近5年內得到重視。AMS是最早集成LAMMPS中reax/c的化學反應分析功能的平臺，同樣僅提供基于分子式的分子數量統計，僅通過提供圖形化的數據提升了分析結果的直觀性；其最近兩年的版本集成了德國亞琛大學D?ntgen等6發布的一個反應分析程序，該程序建立在采用270個原子研究甲烷氧化基元反應的基礎之上，主要進步是可以計算基元反應的反應速率，但難以應用于～10000原子規模的復雜分子體系。MAPS平臺的REAXFF ANALYSIS插件則提供了支持LAMMPS中reax/c的模擬結果表格化顯示的功能，可識別分子、物種、反應以及計算分子壽命7。由此可見，商業化的分子動力學模擬平臺所集成的反應分析程序主要基于開源程序或個別研究者的小程序，在反應分析核心能力的進一步發展上仍然欠缺。最近兩年，在ReaxFF MD模擬的反應網絡自動生成方面取得了進展。華東師范大學朱通等構建的反應網絡自動生成程序ReacNetGenerator采用了隱馬爾科夫鏈方法平滑震蕩的基元反應是一個新的嘗試，可顯著提升生成合理反應網絡的效率，已應用于甲烷燃燒的機理研究和4組分RP-3替代燃料氧化的反應網絡生成8。由于要預先手動標注一部分噪聲反應用于隱馬爾科夫鏈相關矩陣的計算，應用于復雜燃料氧化反應網絡的自動生成還有待更多的工作。上海交通大學吳量和孫淮等9則發展了基于優化的鍵級截斷值計算ReaxFF MD模擬的基元反應速率常數、再結合直接關系圖法對反應進行機理簡化的方法，將其應用于氫氣燃燒的ReaxFF MD模擬結果分析，獲得了合理的簡化動力學模型，展現了從ReaxFF MD模擬基元反應直接獲得動力學性質的潛力。總體而言，ReaxFF MD模擬的反應機理分析方法取得了明顯進展。但與快速發展的ReaxFF MD模擬應用相比，專門用于分析ReaxFF MD模擬結果的方法和程序實現策略仍然滯后，特別是針對大規模復雜體系的ReaxFF MD模擬應用的化學反應分析方法仍然是挑戰。

作者課題組近年來面向國家能源利用相關的重大需求，致力于發展大規模ReaxFF MD模擬的方法，建立了國際首個基于GPU并行ReaxFF MD模擬程序系統GMD-Reax，顯著提升了10000–100000原子規模在單GPU計算節點上模擬計算的效率10；基于化學信息學方法建立了ReaxFF MD反應分析與可視化程序系統VARxMD (Visualization and Analysis of ReaxFF Molecular Dynamics)11,12。相比于國際上已有的ReaxFF MD的反應分析工具，VARxMD在處理大規模分子體系模擬結果方面具有領先優勢和獨特性，其反應分析功能建立在3D化學結構對唯一物種識別、物種反應位點識別、鍵類型識別的基礎之上，可基于相鄰時刻之間的成斷鍵信息自動生成完整的化學反應列表，并進行反應位點的2D和3D結構可視化顯示；再者，基于反應物或產物的化學結構、官能團和反應位點的檢索，可進一步對反應路徑進行分類，并圖形化展示反應路徑的演化13。VARxMD的最新發展是特定反應物和特定生成物之間反應網絡的自動生成與可視化14。VARxMD已經在大規模ReaxFF MD模擬高溫熱解和氧化反應過程、揭示其復雜的反應機理應用中發揮了獨特作用15–18。

在應用于復雜過程模擬的反應機理分析時，目前的VARxMD仍存在一定的局限。首先，目前的VARxMD主要用于分析ReaxFF MD模擬體系所蘊含的化學反應信息，并不包含在經典分子動力學模擬中常用的物理過程參數分析，而一些真實的復雜反應過程是物理和化學變化同時伴隨發生。例如煤在高溫熱解過程中，多個大片段焦油自由基分子會重聚形成更大片段的煤焦前驅體分子。不斷變化的煤焦前驅體分子夾雜在煤熱解體系中，要對其結構特征和時空性質進行表征分析存在困難；難以考察煤熱解反應體系中兩個煤顆粒之間的相互作用；輕焦油從煤顆粒微孔的逸出與微孔的大小、分布的關系也難以獲得。此外，目前VARxMD的反應分析功能主要著眼于模擬體系的全局反應，已有的物種檢索功能和3D化學結構拾取功能可以聚焦于特定的反應物或生成物分子，尚不能聚焦于模擬體系的特定反應區域，但局部區域的反應追蹤與分析對一些特定的過程極為重要。例如含能材料在外部刺激條件下會形成局部升溫并可能演變為爆燃爆轟過程的觸發，考察特定熱點區域的物理和化學變化是研究者極為關注的問題；又如模擬分子篩催化反應體系時，研究人員關心反應物種在分子篩催化劑表面的物理吸附和化學反應細節。然而到目前為止，國內外的ReaxFF MD反應分析程序系統尚未聚焦于復雜體系局部區域的化學反應分析，難以回溯和追蹤特定區域的化學反應。這一重要分析能力的缺失，限制了ReaxFF MD應用于更多復雜過程的反應機理探索。

為此，我們基于在ReaxFF MD模擬反應分析與可視化方面自主研發的VARxMD系統，一方面擴展了物理時空性質演化的分析功能，使得VARxMD可分析ReaxFF MD模擬體系中同時包含的物理和化學變化的演化過程；在此基礎上進一步擴展了VARxMD聚焦于模擬體系3D局部區域內化學反應分析追蹤和物理性質的分析能力。本文概要介紹VARxMD所擴展ReaxFF MD局部區域反應追蹤和物理性質動態分析的建立方法，并通過交互繪制矩形選擇局部區域進行反應追蹤來揭示煤熱解體系中煤顆粒之間孔道的作用、以及通過采用球體局部區域計算區域內徑向分布函數變化進而分析對含能材料熱解過程中所形成的富碳團簇結構進行表征，展示本工作擴展方法的應用及其可能的應用前景。

2 ReaxFF MD模擬體系局部區域反應追蹤及物理性質動態分析方法

從聚焦局部反應區域如表面、反應熱點等角度出發，在VARxMD自動分析模擬體系全局化學反應行為的基礎上，我們為VARxMD程序系統擴展了聚焦所模擬體系中局部反應行為、并追蹤其物理性質和化學反應隨時間的動態演化規律的能力。圖1是本工作擴展的聚焦ReaxFF MD模擬體系局部區域反應分析、物理性質與VARxMD全局區域的關系圖。

圖1 VARxMD系統功能擴展：聚焦ReaxFF MD模擬體系局部區域反應追蹤及物理性質動態分析Fig. 1 Schematic diagram for the extension of VARxMD: tracking local reactions and dynamic evolution of physical properties in a 3D area picked up in a simulation cell of reactive molecular dynamics.

如圖1所示，針對ReaxFF MD模擬體系，VARxMD擴展的局部區域分析模塊建立在已有的全局化學反應分析和物理性質動態分析模塊之上。基于VARxMD全局反應分析和全局的徑向分布函數RDF (Radial Distribution Function)、均方位移MSD (Mean Square Displacement)分析，局部區域分析包括指定局部區域內唯一物種間化學反應的分析追蹤、以及針對特定原子和分子的物理性質分析。

要聚焦局部區域的分析，首先需要對模擬體系中特定的空間進行交互式指定，即進行局部區域空間的選取；接著程序系統自動識別區域內的原子或分子片段；再進行局部區域內以基于物種唯一性的分子片段之間化學反應的識別和分析追蹤、及針對該區域內原子和分子的物理性質分析。局部區域的分析與VARxMD已有的全局分析功能平臺密切相關，且其實現策略直接關系到應用于大規模模擬體系的可行性。局部區域空間的選取本質是ReaxFF MD模擬體系的模擬盒子內一個3D子空間的指定與選取，因此局部區域內的化學反應追蹤與物理時空性質分析建立在VARxMD已有的模擬體系全局3D視圖之上。3D局部區域分析的實現主要包括3個部分：(1)3D局部區域空間的選擇與繪制；(2)局部區域內原子、分子、物種和反應的識別；(3)局部區域化學反應的追蹤與物理性質動態分析。

2.1 ReaxFF MD模擬體系3D局部區域選擇繪制與分子識別

要分析模擬體系中局部空間內復雜物理變化與化學反應隨時間的演化規律，首先要進行3D局部區域的選擇和繪制。在局部區域化學反應追蹤與物理性質動態分析中，針對不同ReaxFF MD模擬可能的應用，如表界面反應、局部的反應熱點、局部的團簇生成反應等，設計和實現了兩類選擇和繪制3D局部區域的模式。一類是用戶在模擬盒子中指定3D局部空間的參考點并交互輸入其幾何參數后自動繪制幾何體的幾何模式，幾何體的形狀目前包括長方體和球體；由于選擇與繪制區域方法是基于面向對象的C++編程模式和模塊化編程實現的，若要擴充新的局部區域選擇幾何體例如圓柱體，可方便地加以實現。另一類交互模式是用戶在3D視圖上利用鼠標繪制任意大小的矩形框，經過光線投射算法識別矩形框后面的實物，再根據實物位置數據繪制長方體。由于VARxMD的3D視圖可視化基于VTK11框架開發，從3D視圖的指定局部區域內識別區域大小并獲取數據的過程可映射為一種消息響應機制。第一類區域選擇方法采用VTK可視化管線機制自動渲染出幾何體空間，其偽代碼如圖2所示。為使VTK窗口交互得到區域數據的渲染功能滿足“高內聚低耦合”的設計要求，第二類局部區域選擇方法采用觀察者設計模式(Observer Design Pattern)與命令設計模式(Command Design Pattern)相結合的策略加以實現，該模式可高效地交互拾取3D視圖中的選中區域對象。

圖2 VARxMD從模擬盒子選取長方體/球體局部區域并識別其中的分子的算法偽代碼Fig. 2 Pseudo code of VARxMD for picking up a 3D area of cuboid/sphere in a simulation cell and for the identification of molecules therein.

VARxMD擴展的第二類局部區域選擇方法與第一類相比，更加靈活。用戶可在3D可視化界面上利用鼠標繪制任意大小的矩形、由程序系統通過光線透射算法識別出射線與矩形框相交的分子。由于用戶在電腦屏幕上選擇的是2D平面，該平面背后的空間實際上并不可見。為此，第二類局部區域選擇的難點在于如何通過光線透射算法選擇出2D平面后的矩形區域、及識別該區域內所有的分子，并對其進行高亮顯示。VTK框架中的觀察者設計模式可讓觀察者相關的目標響應同一個動作指令；命令設計模式可封裝命令，并讓選中的所有目標作為命令接受方，執行VARxMD用戶通過屏幕鼠標觸發的具體命令及操作。結合VTK框架中的觀察者和命令設計模式，VARxMD可高效選中用戶指定的局部區域和該區域內包含的所有分子。其實現策略如圖3所示。

圖3 VARxMD基于觀察者與命令模式結合實現鼠標在屏幕交互繪制矩形選擇3D局部區域實現策略Fig. 3 Implementation strategy for 3D zone pick up interactively in a ReaxFF MD simulation cell using a combination of Observer pattern and Command pattern.

如圖3所示，在VTK 3D可視化編程環境19下，VARxMD通過VTK可視化客戶端模塊監聽和接收真實用戶的屏幕矩形繪制請求。基于觀察者模式，將區域選擇與平面矩形繪制操作建立綁定關系，選擇指定區域，獲得區域幾何體的相關數據；接著利用命令設計模式，接受從該區域背后識別出的分子數據，并對選擇出的分子改變顏色；然后沿著可視化管線的逆方向逐級更新到VTK分子數據源，將VTK分子數據源索引、繼承并搜索由VARxMD全局反應分析獲得的全局分子列表，從而實現了指定局部區域內的分子識別。局部區域分子識別的實現是進行指定局部區域反應分析和追蹤的基礎和關鍵。

2.2 3D局部區域唯一物種間的反應分析與追蹤

當從ReaxFF MD模擬盒子中交互地繪制和選擇了指定的3D局部區域、并對該區域的分子進行識別之后，就可以對該局部區域內當前時刻t的物種(即每一個區分了同分異構體的3D分子結構)與某個時刻t + ?t(當?t> 0為正向追蹤；當?t< 0為逆向回溯)之間的化學反應進行分析和追蹤。圖4給出了VARxMD的3D局部區域反應分析與追蹤的算法偽代碼，其中輸入是從指定3D局部區域中識別出的分子和VARxMD已分析獲得的模擬體系完整反應列表，輸出為當前時刻t至t + ?t之間指定局部區域所有分子涉及的化學反應。如圖4所示，先遍歷?t中每兩個相鄰時間步的時段，在每個時段遍歷當前時刻局部區域所識別出的分子集合，判定分子處于所選擇局部區域的判據是該分子的質心位于該局部區域內，由此得到每個分子及其所屬原子集合；再遍歷原子集合，得到每個原子在相鄰時刻的所屬分子，并對所屬分子去重；最后以當前分子和下一個(或前一個)時刻的分子聯合為索引，在全局反應列表中查找化學反應，如此循環，即可獲得局部區域的化學反應列表。

圖4 VARxMD對ReaxFF MD模擬盒子中局部區域分子的化學反應分析與追蹤算法偽代碼Fig. 4 VARxMD pseudo code of reaction analysis and tracking for a picked 3D area in ReaxFF MD simulation cell.

2.3 3D局部區域物理性質RDF動態分析

VARxMD局部區域物理性質動態分析是指從物理性質的角度分析ReaxFF MD模擬體系中指定局部區域的瞬時結構特征及其演化行為。通過引入經典分子動力學模擬靜態結構的徑向分布函數和均方位移函數等時空性質的計算方法，定量描述所模擬的反應分子體系空間密度隨時間的演化和自擴散行為等。其中，本工作以RDF徑向分布函數為例(局部區域原子MSD的算法請見Supporting Information)，設計并實現了區分原子和分子、并結合多條件同時分析的方法。在指定局部區域的前提下，原子的RDF計算包括不區分元素的所有原子類型、指定元素類型的原子對；并可設置不同時刻以獲得RDF隨時間的動態演化規律。為了方便比較不同元素的原子在特定時刻的結構特征，VARxMD特別設計了同時計算特定時刻、指定不同元素類型的原子RDF和可視化比較的功能。多曲線的RDF繪制變量可以是時間、球殼厚度、某元素原子，能清晰揭示該指定區域內物理結構在時空尺度的演化。

圖5給出了RDF分析設置界面示例，用于分析某含能材料熱分解體系在1250 ps、指定的某個局部區域(區域數據保存為1250 ps. RegAtoms文件)的RDF。此設置將同時分析C-C、C-H、C-O、C-N四種原子對的RDF，球殼間隔為0.02 nm。區域不同元素原子對之間RDF計算可表征該體系指定的局部區域當前時刻的結構特征，分析結果可見3.2節。

圖5 VARxMD對ReaxFF MD模擬的某含材料熱分解體系局部區域在1250 ps分析不同原子對之間RDF計算及多曲線可視化比較的設置界面示例Fig. 5 Screenshot of VARxMD for RDF calculation of varied atom pairs at moment of 1250 ps in a picked 3D area in thermal decomposition of an energetic material system simulated using ReaxFF MD.

3 結果與討論

本文在針對已有的VARxMD全局化學反應分析與可視化系統的基礎上，基于VARxMD中的3D視圖和可視化編程環境VTK，擴展了聚焦ReaxFF MD模擬體系中指定局部區域的反應追蹤及物理性質動態分析功能。此擴展為ReaxFF MD模擬應用于復雜反應過程、揭示重點區域的化學細節和探測其物理性質提供了可能。下面兩節概述該新功能應用于研究煤熱解過程中煤顆粒間孔道的反應分析和含能材料熱分解過程中所形成納米碳團簇的物理性質。

3.1 VARxMD的3D局部區域反應追蹤應用于煤熱解過程煤顆粒孔道的相互作用

煤熱解是煤的熱加工利用如氣化、燃燒的基礎過程，了解煤熱解機理可為煤清潔利用提供理論支持。由于煤熱解是自由基驅動的高溫快速反應過程，已有的實驗手段難以原位捕捉煤熱解過程中的中間體和反應細節。基于對大規模煤分子模型的ReaxFF MD模擬可揭示煤熱解過程反應、中間體及熱解產物的演化規律，是一種有潛力深入認識煤熱解本征反應的方法20–22。

一般認為，煤熱解過程起始于煤結構中最弱的橋鍵，橋鍵斷裂之后，某些較小的揮發分中間體會遷移至游離空間發生二次裂解或者被小分子自由基(HO、CH3等)穩定下來，促進煤熱解過程23。因此，兩個煤顆粒之間的孔道(游離空間)發生的化學反應對認識煤熱解過程至關重要。由于反應分析能力的限制，已有的煤熱解反應模擬主要針對單個煤顆粒內部的模擬。3D局部區域反應分析與追蹤為研究兩個煤顆粒孔道之間的反應細節提供了可能。圖6a是利用VARxMD新功能分析的一個包含了兩個煤顆粒的大規模煤模型在高溫的熱解反應的示例。該煤模型含99544個原子，分子式為C46920H43320O9088N72S144，其中兩煤顆粒間的孔道間距約為1 nm。利用ReaxFF MD，采用慢速升溫速率2 K?ps?1，將該煤模型從500 K加熱至2500 K。本次示例以10 ps為時間間隔，分析了從5–205 ps(510–910 K)、共20幀模擬軌跡的初始煤熱解過程，共得到4782個化學反應。

圖6 (a)VARxMD的屏幕交互繪制矩形選擇3D局部區域應用于ReaxFF MD模擬的煤熱解體系煤顆粒孔道間化學反應追蹤分析的屏幕截圖示例以及(b)5 ps和205 ps 3D視圖的對比Fig. 6 (a)VARxMD screenshot of the 3D picked zone and the reactions tracking therein through interactive drawing of a rectangle on screen, focused on channel reactions and (b)comparison of 3D views at 5 ps and 205 ps of coal particle thermal systems simulated using ReaxFF MD.

如圖6所示，利用VARxMD新擴展的屏幕交互繪制矩形模式(第二類選擇模式)選擇兩個煤顆粒之間的孔道空間，即圖中3D窗口紅色邊框的長方體局部區域。該長方體沿X軸方向的長度為2.612 nm，占模擬盒子X軸向長度20 nm的13.06%。長方體內藍色高亮部分為選中的3D局部區域內所識別出的原子與鍵；右側是以樹圖結構展示的局部區域原子所參與的化學反應列表，每個反應的細節可在反應列表樹圖下方以2D化學結構式清晰展示，也可通過Region Species Generator檢索識別出該局部區域內所有物種的分子集合，并基于元素、分子式、基團等進行反應物或生成物的統計和歸類。圖6b給出了指定局部區域在5 ps和205 ps的3D視圖對比，可清晰看出隨著煤顆粒間的反應發生，局部區域逐漸被熱解反應生成的產物和中間體分子占據，煤顆粒間的較大孔道因劇烈反應發生變小，微孔道增多。在5–205 ps可分析得到局部區域的化學反應為932個。

表1給出了在5–205 ps的煤熱解初始過程中指定局部區域發生的化學反應數量與全局反應數量的對比。如表1所示，雖然所選取的孔道局部區域體積只占整個煤模型體系體積的13.06%，該區域發生的化學反應在裂解起始階段占了全局體系的20%左右，說明該區域的確是煤熱解的活躍區域，生成的中間體和自由基會促進整個煤分子的熱解，增加煤焦油的收率。通過本工作實現的局部區域反應細節追蹤的獨特功能，可獲得在該段模擬時間獲得的特征化學反應及其2D結構。表2舉例列出了部分反應的化學方程式，基于2D結構的化學反應表示在Supporting Information中提供，這些反應物種的2D結構信息可幫助深入認識該區域的相互作用和所發生的化學變化。如表2所示，該區域發生了大部分煤熱解的特征反應，包括最薄弱的－O－CH2－橋鍵斷裂(表2中的反應2、5、6、7、9、12、16、17)，HO、CH3等游離小分子自由基的生成(表2中的反應6、8、14)，煤活化的結構轉化和可逆震蕩反應(表2中的反應3、4)，以及由于HO、CHO2脫落和H轉移反應促進的H2O和CO2生成(表2中的反應10、15、16、17)。此外，由于孔道空間較大，有利于小分子自由基的移動，該局部區域也發生了單一煤顆粒分子模擬時沒有發現的初始裂解過程中H2的生成(表2中的反應1)和輕質焦油的生成(表2中的反應11)。這些局部區域的反應細節有利于幫助認識孔道對煤熱解的作用，是其它分析軟件或擴展之前的VARxMD功能無法提供的。

表1 圖6煤熱解模擬體系內煤顆粒孔道局部區域化學反應與全局模擬體系反應在5–205 ps間的數量對比Table 1 Number comparison of the local chemical reactions in the picked zone with the total reactions in thecoal pyrolysis simulation system of Fig. 6 during the period of 5–205 ps.

表2 5–205 ps模擬時間內煤顆粒體系特征化學反應列表Table 2 Examples of characteristic reaction of coal particle system during 5–205 ps.

3.2 VARxMD的3D局部區域物理性質動態分析應用于含能材料熱解過程的團簇生成

近年來，基于ReaxFF力場的反應分子動力學方法在認識含能材料復雜反應行為的研究非常活躍，從微觀尺度揭示了多種含能材料在熱載荷24、壓縮25、沖擊26,27以及剪切28等不同外部刺激下的化學響應機制，而且為了解炸藥在高溫高壓極端爆轟條件下復雜化學過程提供了可能。團簇大分子是富碳含能材料TNT爆轟過程中的重要產物，與炸藥的爆轟過程和做功能力密切相關，是含能材料爆轟領域的研究熱點29。以利用ReaxFF MD模擬研究大規模TNT晶體模型(經驗化學式為C8960H6400O7680N3840)在高溫3000 K的熱分解-膨脹-冷卻過程為例，借助VARxMD最新擴展的局部區域物理性質動態分析功能，可幫助認識TNT高溫熱分解過程中富碳團簇的生成機理。該模擬體系含1280個TNT分子，共26880個原子，模擬的總時長為1350 ps，時間步長0.1 fs，軌跡輸出間隔為100 fs。利用VARxMD以較大的時間間隔25 ps分析完整的模擬軌跡，可獲得5150個反應。

在TNT高溫熱分解模擬的末期可觀察到許多大小不同的原子團簇形成(如圖7左側所示，模擬盒子在1250 ps時刻的3D截圖)。為細致觀察所生成原子團簇，在VARxMD的3D可視化窗口中交互選擇了模擬盒子中一處包含層狀結構的球形區域(半徑為2 nm，如圖7右側所示)作為重點分析的對象。該局部區域內包含幾層大富碳分子碎片及一些CO2，N2，H2O，H2等小分子氣體產物。為進一步分析該層狀團簇的結構特征，計算1250 ps時刻下該局部球形區域內不同原子的徑向分布函數。首先分析了選中局部球形區域內C-C原子對的RDF曲線(間隔參數為0.002 nm)以了解原子團簇碳環骨架的結構特征，如圖8a所示。與5層石墨烯結構的C-C原子對RDF進行對比，可以看到該層狀分子團簇的C-C原子對的RDF主要峰值位置和石墨烯吻合較好，這表明該層狀分子團簇具有類石墨烯結構。觀察該層狀原子團簇各片層分子的結構(如圖8b所示)，可以看到每一分子碎片的主體骨架均具有稠環結構，以六元碳環為主、夾雜少量五元和七元碳環結構，每一分子片層的周圍還帶有一些側鏈基團。這些結構特點說明該層狀分子團簇尚未形成完美的石墨烯六元稠環結構，仍需更長時間的演化。

圖7 VARxMD的繪制球形選擇3D局部區域功能應用于ReaxFF MD模擬TNT熱解體系在1250 ps的界面截圖Fig. 7 VARxMD interface screenshot of drawing a sphere and selecting 3D local area function applied to ReaxFF MD simulation TNT pyrolysis system at 1250 ps.

圖8 (a)TNT晶體模型熱解體系的球形區域層狀分子團C-C的RDF與5層石墨烯C-C的標準RDF對比及(b)層狀團簇中主要大分子結構Fig. 8 (a)Comparison of the RDF of the spherical region layered molecular clusters C-C of the pyrolysis system of the TNT crystal model with the standard RDF of the 5-layer graphene standard C-C and(b)major macromolecular structures in layered clusters.

選中的層狀結構的大分子團簇中還含有一些N、O雜原子，這與已報道的富碳含能材料爆轟得到碳納米材料中常含有N、O雜原子是吻合的30。因此利用新擴展的局部區域多條件RDF分析功能，計算了區域內其他原子對C-O、C-N和C-H之間的RDF(間隔參數為0.02 nm)，如圖9所示。由不同原子的徑向分布函數曲線的最大峰值可以看到，C－O、C－N和C－H鍵的鍵長分別在0.13 nm、0.13 nm和0.11 nm左右。與標準原子間共價鍵鍵長比較，可確定團簇中C和O原子之間多為單鍵，而C和N原子間主要為雙鍵。

圖9 VARxMD分析ReaxFF MD模擬的TNT晶體模型熱解體系的球形局部區域在1250 ps，0.02 nm下的不同原子對之間RDF計算結果的顯示界面截圖Fig. 9 VARxMD analysis screenshot of the display interface of the RDF calculation results between different atomic pairs at 1250 ps, 0.02 nm in the spherical partial region of the TNT crystal model pyrolysis system simulated by ReaxFF MD.

4 結論

本工作基于自主研發的、先進的反應分子動力學模擬結果反應分析與可視化工具VARxMD程序系統，以VARxMD的3D可視化視圖作為切入點，為VARxMD進一步擴展了對指定的3D局部區域內的進行反應追蹤與瞬態物理性質分析能力。從ReaxFF MD模擬體系全局分析到3D局部區域分析能力的擴展使得研究者聚焦分析模擬體系內、特定的局部反應熱點成為可能，已應用于煤熱解反應模擬中煤顆粒間的反應追蹤以及含能材料熱分解過程中富碳團簇形成中的瞬態結構表征，也有望應用于催化反應體系表界面反應分析、含能材料爆轟過程中反應熱點分析等復雜反應過程的研究。

Supporting Information:available free of chargeviathe internet at http://www.whxb.pku.edu.cn.