ReaxFF MD局部區(qū)域反應(yīng)追蹤與物理性質(zhì)可視化分析

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

1中國科學(xué)院過程工程研究所多相復(fù)雜系統(tǒng)國家重點(diǎn)實(shí)驗(yàn)室，北京 100190

2中國科學(xué)院大學(xué)化學(xué)工程學(xué)院，北京 100049

3中國科學(xué)院綠色制造創(chuàng)新研究院，北京 100190

關(guān)鍵字：ReaxFF MD；區(qū)域反應(yīng)追蹤；物理性質(zhì)；可視化

1 引言

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

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

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

作者課題組近年來面向國家能源利用相關(guān)的重大需求，致力于發(fā)展大規(guī)模ReaxFF MD模擬的方法，建立了國際首個(gè)基于GPU并行ReaxFF MD模擬程序系統(tǒng)GMD-Reax，顯著提升了10000–100000原子規(guī)模在單GPU計(jì)算節(jié)點(diǎn)上模擬計(jì)算的效率10；基于化學(xué)信息學(xué)方法建立了ReaxFF MD反應(yīng)分析與可視化程序系統(tǒng)VARxMD (Visualization and Analysis of ReaxFF Molecular Dynamics)11,12。相比于國際上已有的ReaxFF MD的反應(yīng)分析工具，VARxMD在處理大規(guī)模分子體系模擬結(jié)果方面具有領(lǐng)先優(yōu)勢(shì)和獨(dú)特性，其反應(yīng)分析功能建立在3D化學(xué)結(jié)構(gòu)對(duì)唯一物種識(shí)別、物種反應(yīng)位點(diǎn)識(shí)別、鍵類型識(shí)別的基礎(chǔ)之上，可基于相鄰時(shí)刻之間的成斷鍵信息自動(dòng)生成完整的化學(xué)反應(yīng)列表，并進(jìn)行反應(yīng)位點(diǎn)的2D和3D結(jié)構(gòu)可視化顯示；再者，基于反應(yīng)物或產(chǎn)物的化學(xué)結(jié)構(gòu)、官能團(tuán)和反應(yīng)位點(diǎn)的檢索，可進(jìn)一步對(duì)反應(yīng)路徑進(jìn)行分類，并圖形化展示反應(yīng)路徑的演化13。VARxMD的最新發(fā)展是特定反應(yīng)物和特定生成物之間反應(yīng)網(wǎng)絡(luò)的自動(dòng)生成與可視化14。VARxMD已經(jīng)在大規(guī)模ReaxFF MD模擬高溫?zé)峤夂脱趸磻?yīng)過程、揭示其復(fù)雜的反應(yīng)機(jī)理應(yīng)用中發(fā)揮了獨(dú)特作用15–18。

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

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

2 ReaxFF MD模擬體系局部區(qū)域反應(yīng)追蹤及物理性質(zhì)動(dòng)態(tài)分析方法

從聚焦局部反應(yīng)區(qū)域如表面、反應(yīng)熱點(diǎn)等角度出發(fā)，在VARxMD自動(dòng)分析模擬體系全局化學(xué)反應(yīng)行為的基礎(chǔ)上，我們?yōu)閂ARxMD程序系統(tǒng)擴(kuò)展了聚焦所模擬體系中局部反應(yīng)行為、并追蹤其物理性質(zhì)和化學(xué)反應(yīng)隨時(shí)間的動(dòng)態(tài)演化規(guī)律的能力。圖1是本工作擴(kuò)展的聚焦ReaxFF MD模擬體系局部區(qū)域反應(yīng)分析、物理性質(zhì)與VARxMD全局區(qū)域的關(guān)系圖。

圖1 VARxMD系統(tǒng)功能擴(kuò)展：聚焦ReaxFF MD模擬體系局部區(qū)域反應(yīng)追蹤及物理性質(zhì)動(dòng)態(tài)分析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所示，針對(duì)ReaxFF MD模擬體系，VARxMD擴(kuò)展的局部區(qū)域分析模塊建立在已有的全局化學(xué)反應(yīng)分析和物理性質(zhì)動(dòng)態(tài)分析模塊之上?；赩ARxMD全局反應(yīng)分析和全局的徑向分布函數(shù)RDF (Radial Distribution Function)、均方位移MSD (Mean Square Displacement)分析，局部區(qū)域分析包括指定局部區(qū)域內(nèi)唯一物種間化學(xué)反應(yīng)的分析追蹤、以及針對(duì)特定原子和分子的物理性質(zhì)分析。

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

2.1 ReaxFF MD模擬體系3D局部區(qū)域選擇繪制與分子識(shí)別

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

圖2 VARxMD從模擬盒子選取長方體/球體局部區(qū)域并識(shí)別其中的分子的算法偽代碼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擴(kuò)展的第二類局部區(qū)域選擇方法與第一類相比，更加靈活。用戶可在3D可視化界面上利用鼠標(biāo)繪制任意大小的矩形、由程序系統(tǒng)通過光線透射算法識(shí)別出射線與矩形框相交的分子。由于用戶在電腦屏幕上選擇的是2D平面，該平面背后的空間實(shí)際上并不可見。為此，第二類局部區(qū)域選擇的難點(diǎn)在于如何通過光線透射算法選擇出2D平面后的矩形區(qū)域、及識(shí)別該區(qū)域內(nèi)所有的分子，并對(duì)其進(jìn)行高亮顯示。VTK框架中的觀察者設(shè)計(jì)模式可讓觀察者相關(guān)的目標(biāo)響應(yīng)同一個(gè)動(dòng)作指令；命令設(shè)計(jì)模式可封裝命令，并讓選中的所有目標(biāo)作為命令接受方，執(zhí)行VARxMD用戶通過屏幕鼠標(biāo)觸發(fā)的具體命令及操作。結(jié)合VTK框架中的觀察者和命令設(shè)計(jì)模式，VARxMD可高效選中用戶指定的局部區(qū)域和該區(qū)域內(nèi)包含的所有分子。其實(shí)現(xiàn)策略如圖3所示。

圖3 VARxMD基于觀察者與命令模式結(jié)合實(shí)現(xiàn)鼠標(biāo)在屏幕交互繪制矩形選擇3D局部區(qū)域?qū)崿F(xiàn)策略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可視化編程環(huán)境19下，VARxMD通過VTK可視化客戶端模塊監(jiān)聽和接收真實(shí)用戶的屏幕矩形繪制請(qǐng)求?；谟^察者模式，將區(qū)域選擇與平面矩形繪制操作建立綁定關(guān)系，選擇指定區(qū)域，獲得區(qū)域幾何體的相關(guān)數(shù)據(jù)；接著利用命令設(shè)計(jì)模式，接受從該區(qū)域背后識(shí)別出的分子數(shù)據(jù)，并對(duì)選擇出的分子改變顏色；然后沿著可視化管線的逆方向逐級(jí)更新到VTK分子數(shù)據(jù)源，將VTK分子數(shù)據(jù)源索引、繼承并搜索由VARxMD全局反應(yīng)分析獲得的全局分子列表，從而實(shí)現(xiàn)了指定局部區(qū)域內(nèi)的分子識(shí)別。局部區(qū)域分子識(shí)別的實(shí)現(xiàn)是進(jìn)行指定局部區(qū)域反應(yīng)分析和追蹤的基礎(chǔ)和關(guān)鍵。

2.2 3D局部區(qū)域唯一物種間的反應(yīng)分析與追蹤

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

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

2.3 3D局部區(qū)域物理性質(zhì)RDF動(dòng)態(tài)分析

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

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

圖5 VARxMD對(duì)ReaxFF MD模擬的某含材料熱分解體系局部區(qū)域在1250 ps分析不同原子對(duì)之間RDF計(jì)算及多曲線可視化比較的設(shè)置界面示例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 結(jié)果與討論

本文在針對(duì)已有的VARxMD全局化學(xué)反應(yīng)分析與可視化系統(tǒng)的基礎(chǔ)上，基于VARxMD中的3D視圖和可視化編程環(huán)境VTK，擴(kuò)展了聚焦ReaxFF MD模擬體系中指定局部區(qū)域的反應(yīng)追蹤及物理性質(zhì)動(dòng)態(tài)分析功能。此擴(kuò)展為ReaxFF MD模擬應(yīng)用于復(fù)雜反應(yīng)過程、揭示重點(diǎn)區(qū)域的化學(xué)細(xì)節(jié)和探測(cè)其物理性質(zhì)提供了可能。下面兩節(jié)概述該新功能應(yīng)用于研究煤熱解過程中煤顆粒間孔道的反應(yīng)分析和含能材料熱分解過程中所形成納米碳團(tuán)簇的物理性質(zhì)。

3.1 VARxMD的3D局部區(qū)域反應(yīng)追蹤應(yīng)用于煤熱解過程煤顆?？椎赖南嗷プ饔?/h3>

煤熱解是煤的熱加工利用如氣化、燃燒的基礎(chǔ)過程，了解煤熱解機(jī)理可為煤清潔利用提供理論支持。由于煤熱解是自由基驅(qū)動(dòng)的高溫快速反應(yīng)過程，已有的實(shí)驗(yàn)手段難以原位捕捉煤熱解過程中的中間體和反應(yīng)細(xì)節(jié)?；趯?duì)大規(guī)模煤分子模型的ReaxFF MD模擬可揭示煤熱解過程反應(yīng)、中間體及熱解產(chǎn)物的演化規(guī)律，是一種有潛力深入認(rèn)識(shí)煤熱解本征反應(yīng)的方法20–22。

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

圖6 (a)VARxMD的屏幕交互繪制矩形選擇3D局部區(qū)域應(yīng)用于ReaxFF MD模擬的煤熱解體系煤顆?？椎篱g化學(xué)反應(yīng)追蹤分析的屏幕截圖示例以及(b)5 ps和205 ps 3D視圖的對(duì)比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新擴(kuò)展的屏幕交互繪制矩形模式(第二類選擇模式)選擇兩個(gè)煤顆粒之間的孔道空間，即圖中3D窗口紅色邊框的長方體局部區(qū)域。該長方體沿X軸方向的長度為2.612 nm，占模擬盒子X軸向長度20 nm的13.06%。長方體內(nèi)藍(lán)色高亮部分為選中的3D局部區(qū)域內(nèi)所識(shí)別出的原子與鍵；右側(cè)是以樹圖結(jié)構(gòu)展示的局部區(qū)域原子所參與的化學(xué)反應(yīng)列表，每個(gè)反應(yīng)的細(xì)節(jié)可在反應(yīng)列表樹圖下方以2D化學(xué)結(jié)構(gòu)式清晰展示，也可通過Region Species Generator檢索識(shí)別出該局部區(qū)域內(nèi)所有物種的分子集合，并基于元素、分子式、基團(tuán)等進(jìn)行反應(yīng)物或生成物的統(tǒng)計(jì)和歸類。圖6b給出了指定局部區(qū)域在5 ps和205 ps的3D視圖對(duì)比，可清晰看出隨著煤顆粒間的反應(yīng)發(fā)生，局部區(qū)域逐漸被熱解反應(yīng)生成的產(chǎn)物和中間體分子占據(jù)，煤顆粒間的較大孔道因劇烈反應(yīng)發(fā)生變小，微孔道增多。在5–205 ps可分析得到局部區(qū)域的化學(xué)反應(yīng)為932個(gè)。

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

表1 圖6煤熱解模擬體系內(nèi)煤顆?？椎谰植繀^(qū)域化學(xué)反應(yīng)與全局模擬體系反應(yīng)在5–205 ps間的數(shù)量對(duì)比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模擬時(shí)間內(nèi)煤顆粒體系特征化學(xué)反應(yīng)列表Table 2 Examples of characteristic reaction of coal particle system during 5–205 ps.

3.2 VARxMD的3D局部區(qū)域物理性質(zhì)動(dòng)態(tài)分析應(yīng)用于含能材料熱解過程的團(tuán)簇生成

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

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

圖7 VARxMD的繪制球形選擇3D局部區(qū)域功能應(yīng)用于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晶體模型熱解體系的球形區(qū)域?qū)訝罘肿訄F(tuán)C-C的RDF與5層石墨烯C-C的標(biāo)準(zhǔn)RDF對(duì)比及(b)層狀團(tuán)簇中主要大分子結(jié)構(gòu)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.

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

圖9 VARxMD分析ReaxFF MD模擬的TNT晶體模型熱解體系的球形局部區(qū)域在1250 ps，0.02 nm下的不同原子對(duì)之間RDF計(jì)算結(jié)果的顯示界面截圖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 結(jié)論

本工作基于自主研發(fā)的、先進(jìn)的反應(yīng)分子動(dòng)力學(xué)模擬結(jié)果反應(yīng)分析與可視化工具VARxMD程序系統(tǒng)，以VARxMD的3D可視化視圖作為切入點(diǎn)，為VARxMD進(jìn)一步擴(kuò)展了對(duì)指定的3D局部區(qū)域內(nèi)的進(jìn)行反應(yīng)追蹤與瞬態(tài)物理性質(zhì)分析能力。從ReaxFF MD模擬體系全局分析到3D局部區(qū)域分析能力的擴(kuò)展使得研究者聚焦分析模擬體系內(nèi)、特定的局部反應(yīng)熱點(diǎn)成為可能，已應(yīng)用于煤熱解反應(yīng)模擬中煤顆粒間的反應(yīng)追蹤以及含能材料熱分解過程中富碳團(tuán)簇形成中的瞬態(tài)結(jié)構(gòu)表征，也有望應(yīng)用于催化反應(yīng)體系表界面反應(yīng)分析、含能材料爆轟過程中反應(yīng)熱點(diǎn)分析等復(fù)雜反應(yīng)過程的研究。

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