劉英哲, 來蔚鵬, 王 玉, 尉 濤, 任 淦, 葛忠學, 康 瑩

(西安近代化學研究所氟氮化工資源高效開發(fā)與利用國家重點實驗室, 陜西 西安 710065)

1 引 言

石墨烯(GRA)具有優(yōu)異的結(jié)構(gòu)、物化及力學特性,在含能材料領(lǐng)域具有廣闊的應用前景,如降低機械感度,提高熱分解穩(wěn)定性,增強能量釋放效率等[1]。曾有研究表明，功能化的GRA對硝基甲烷(NM)的燃燒反應具有催化作用,能夠提高NM的燃燒速率[2]。隨后,基于第一性原理與反應力場的分子動力學模擬分別揭示出該催化效應的機理源于NM與GRA表面官能團發(fā)生的原子交換[3-4]。隨著反應溫度的增加,GRA也顯示出了與功能化GRA相似的催化活性[4]。另外,量子化學計算結(jié)果顯示GRA表面能夠影響NM三種初始反應活化能的大小,并改變了活化能的順序,也從一定程度上體現(xiàn)了GRA的催化作用[5]。然而,關(guān)于GRA對NM結(jié)構(gòu)排布的影響研究卻較少。

密度泛函理論研究結(jié)果表明NM分子間相對取向的改變可以引起CN鍵離解能的變化,導致NM在不同的外界條件刺激下呈現(xiàn)出不同的敏感程度[6],這種由結(jié)構(gòu)排布差異所引起的性能變化(如感度、穩(wěn)定性等)在炸藥晶體中表現(xiàn)尤為明顯。因此,在GRA平面結(jié)構(gòu)的誘導下,NM的結(jié)構(gòu)排布很可能與液體狀態(tài)下的特征不同,導致分子間的反應傳遞過程發(fā)生變化,從而影響最終的反應性能。為此,本研究構(gòu)建了NM-GRA模型,借助分子動力學模擬方法從分子尺度上模擬了NM在GRA表面的動力學行為,通過分析密度分布、偶極矩取向、有序性等,揭示了NM在GRA表面的結(jié)構(gòu)排布特征,旨在為基于結(jié)構(gòu)排布差異理解GRA對NM燃燒的催化效應提供一個基本前提。

2 計算方法

2.1 結(jié)構(gòu)模型

NM溶劑的初始結(jié)構(gòu)取自先前的研究工作[7]。采用CHARMM36分子力場描述NM分子,該力場能較好地重現(xiàn)出液體NM的密度[8]。通過VMD1.9.1軟件[9]中Carbon Nanostructure Builder插件構(gòu)建單層理想的扶手椅GRA結(jié)構(gòu),平面尺寸約為25 ?×25 ?,并將其法向量沿z軸擺放(見圖1)。采用CHARMM27分子力場[10]中的sp2碳原子參數(shù)描述GRA,原子點電荷設置為0 e。整個模擬體系共包含10351個原子,對應的模擬盒子尺寸約為51 ?×51 ?×51 ?。為了模擬與數(shù)據(jù)分析的方便,借助1.0 kcal/mol/?2的弱簡諧勢將GRA限制在初始位置附近。

圖1 NM與GRA模擬體系的初始結(jié)構(gòu)

Fig.1 Initial structure of a GRA sheet immersed in liquid NM solvent for the MD simulation

2.2 分子動力學模擬

采用軟件NAMD2.9[11]在NPT(等溫等壓)系綜中進行分子動力學模擬,在笛卡爾坐標系的三維方向上施加周期性邊界條件。采用SHAKE/RATTLE方法[12-13]將含有氫原子的共價鍵長度限制在平衡值。范德華作用的截斷半徑取12 ?,長程靜電作用由PME算法[14]計算。通過多步Verlet r-RESPA算法[15]對牛頓方程進行積分,短程作用的積分步長設置為2 fs,長程作用的積分步長設置為4 fs。借助Langevin動力學和Langevin控壓方法[16]將溫度和壓強分別控制在300 K和105Pa。體系在最終模擬前經(jīng)歷了如下的平衡步驟: (1) 5000步的能量最小化; (2) 從0 K到300 K的50 ps升溫過程; (3) 300 K和1 bar條件下50 ps的平衡過程。然后,將模擬繼續(xù)運行10 ns。軌跡分析及可視化由VMD1.9.1軟件[9]處理。

3 結(jié)果與討論

3.1 相互作用能

為了考察體系在10 ns的模擬時間尺度內(nèi)是否達到充分平衡,監(jiān)測了GRA與NM相互作用能(Einter)隨模擬時間的演化情況，其計算公式如式(1)：

Einter=Etotal-EGRA-ENM

(1)

式中,Etotal,EGRA,ENM分別為整個體系、GRA、NM溶劑的勢能,kJ·mol-1。GRA與NM相互作用能隨模擬時間的演化圖如圖2所示，由圖2可知，Einter在整個模擬過程中平穩(wěn)地波動,說明10 ns的模擬時間尺度對于NM在GRA表面的結(jié)構(gòu)重排過程是足夠的。因此,選擇最后5 ns的模擬軌跡進行后續(xù)的數(shù)據(jù)分析。

圖2 GRA與NM相互作用能隨模擬時間的演化圖

Fig.2 Time evolution of interaction energies between GRA and NM

3.2 密度分布與結(jié)構(gòu)

為了考察GRA表面對NM結(jié)構(gòu)排布的影響,首先統(tǒng)計了NM密度沿GRA平面法向量方向的變化情況,結(jié)果如圖3所示。從圖3可以看出,隨著NM與GRA平面垂直距離的增加,其密度呈波動式衰減,最終趨于液體NM的平均密度值。根據(jù)峰的位置,可以將NM沿GRA平面法方向量方向的分布分為4個區(qū)域,分別記為L1(2.0～5.9 ?)、L2(5.9～9.8 ?)、L3(9.8～13.7 ?)和L4(13.7～17.6 ?),這也說明在GRA平面結(jié)構(gòu)的誘導下,NM呈層狀結(jié)構(gòu)分布。其中,L1層內(nèi)的NM密度最高,接近3.0 g·cm-3,約為液體NM平均密度的2.5倍,說明在GRA范德華作用的驅(qū)動下,大量NM分子聚集在GRA表面。當NM與GRA平面距離增加時,范德華作用隨之減弱,導致NM密度降低,從L3層開始,NM密度已逐漸接近平均密度。

圖3 NM及其官能團沿GRA平面法向量方向的密度分布圖,黑色虛線為液體NM的平均密度

Fig.3 Mass density profiles of NM, methyl group, and nitro group determined along the normal to the GRA surface as a function of distance. The dash line is the average mass density of bulk NM

另外，計算了NM中—NO2與—CH3基團的密度分布，結(jié)果如圖3所示，由圖3可知，在L1層中，—NO2基團存在兩個分布峰，而—CH3基團存在一個分布峰，說明NM分子中—NO2基團相對于—CH3基團存在兩種典型的結(jié)構(gòu)取向。由于—NO2和—CH3基團可以繞C—N鍵進行旋轉(zhuǎn)，導致NM有兩種不同的異構(gòu)體：交錯式與重疊式(見圖4)，對應的H—C—N—O二面角的最小值分別為30°和0°。在真空中，交錯式與重疊式的轉(zhuǎn)化能壘小于0.05 kJ·mol-1，表明—NO2與—CH3基團幾乎能夠繞著C—N鍵自由地旋轉(zhuǎn)[17]。在液體NM中，分子動力學模擬結(jié)果表明兩種結(jié)構(gòu)具有同等的分布概率[7]。然而，當NM在GRA表面時，量化計算顯示NM基本保持交錯式的結(jié)構(gòu)，此時NM能夠最大程度地與GRA平面相互作用[5]。為此，統(tǒng)計了不同層中兩種NM構(gòu)型出現(xiàn)的概率，這里將H—C—N—O二面角的最小值處于0～15°及15～30°的NM結(jié)構(gòu)分別定義為重疊式與交錯式。圖4c為不同層中交錯式與重疊式結(jié)構(gòu)的概率比。從圖4c可以看出，在L1層內(nèi)，交錯式結(jié)構(gòu)出現(xiàn)的概率約為重疊式結(jié)構(gòu)的2倍，而在其他層中，二者出現(xiàn)的概率基本相等，這與圖3中—NO2與—CH3基團的密度分布是一致的，說明在GRA平面結(jié)構(gòu)的誘導下，部分重疊式結(jié)構(gòu)轉(zhuǎn)化為交錯式結(jié)構(gòu)。

a. staggered NM structure b. eclipsed NM structure

c. probability ratio of staggered NM to eclipsed NM in different solvent layers

圖4 NM兩種構(gòu)型及其相對比例

Fig.4 Two structures of MN and their relative ratio

3.3 偶極矩取向

由圖3可知，在GRA表面NM具有很高的局域密度，在如此密集的結(jié)構(gòu)排布中，NM很可能具有某種特定的取向。為了考察NM的排布取向，這里將NM偶極矩矢量與z軸方向的夾角定義為取向角θ(°)。圖5為不同層內(nèi)NM的cosθ分布圖。從圖5中可以看出，在L1層內(nèi)，取向角在cosθ為0處具有明顯的單峰分布，表明NM具有特定的排布取向，即偶極矩矢量平行于GRA平面。在L2層內(nèi)，NM也具有該特定取向，但分布概率遠低于L1層。對于L3和L4層，cosθ分布曲線幾乎為一條直線，說明NM已不再具有特定的排布取向，這與液體NM內(nèi)的分布情況是一致的[7]，同時也表明GRA對L3和L4層內(nèi)NM的結(jié)構(gòu)排布基本沒有影響。

圖5 不同層內(nèi)NM的cosθ分布圖

Fig.5 Distributions of cosθfor different layers

不同層內(nèi)NM偶極矩大小的分布概率如圖6所示。由圖6可知，所有層內(nèi)的偶極矩分布情況幾乎完全一樣，說明GRA平面主要影響了NM的偶極矩取向，并沒有改變偶極矩的大小，這可能與GRA中碳原子點電荷為0有關(guān)，導致GRA與NM之間不存在靜電相互作用，從而難以影響偶極矩的大小。

圖6 不同層內(nèi)NM偶極矩大小的分布圖

Fig.6 Distributions of the NM dipole moments for different layers

3.4 有序性

圖5也體現(xiàn)了NM結(jié)構(gòu)排布的有序程度，即cosθ分布越集中，則對應的結(jié)構(gòu)排布越有序。因此，采用序參數(shù)(S)進一步定量地表示NM結(jié)構(gòu)排布的有序程度[18]，其計算公式如式(2)：

(2)

由式(2)可知，當序參數(shù)分別為1.0，0.0和-0.5時，NM偶極矩矢量平行、隨機及垂直于z軸。圖7為不同層內(nèi)的平均序參數(shù)。從圖7可以看出，L1層的序參數(shù)接近-0.5，結(jié)構(gòu)排布的有序程度最高，說明大部分的NM偶極矩矢量垂直于z軸，即平行于GRA表面。相比之下，NM在L2層內(nèi)的排布也具有一定的有序性，但有序化程度不高。至于L3和L4層，NM的結(jié)構(gòu)排布則不存在有序性，偶極矩矢量的取向完全隨機化。因此，GRA表面能夠誘導NM進行有序化的結(jié)構(gòu)排布，但隨著與GRA表面距離的增加，有序程度逐漸降低。

圖7 不同層內(nèi)NM結(jié)構(gòu)排布的平均序參數(shù)

Fig.7 Average order parameterSfor different layers

4 結(jié) 論

(1) NM在GRA表面呈層狀分布，在L1層內(nèi)，交錯式構(gòu)型出現(xiàn)的概率約為重疊式構(gòu)型的2倍，局域密度大于平均密度，結(jié)構(gòu)排布較為密集。

(2) NM在GRA表面呈有序化排布，L1層的有序性最高，NM具有特定的偶極矩取向，即偶極矩矢量平行于GRA平面，但偶極矩的大小并沒有改變。

(3) 隨著NM與GRA距離的增加，GRA對NM結(jié)構(gòu)排布的影響減弱，L3和L4層基本與液體NM的結(jié)構(gòu)排布一致。

參考文獻:

[1] 蘭元飛, 李霄羽, 羅運軍. 石墨烯在含能材料中的應用研究進展[J]. 火炸藥學報, 2015, 38(1): 1-7.

LAN Yuan-fei, LI Xiao-yu, LUO Yun-jun. Research progress on application of grapheme in energetic materials[J].ChineseJournalofExplosive&Propellants, 2015, 38(1): 1-7.

[2] Sabourin J L, Dabbs D M, Yetter R A, et al. Functionalized graphene sheet colloids for enhanced fuel/propellant combustion[J].ACSNano, 2009, 3(12): 3945-3954.

[3] Liu L M, Car R, Selloni, A, et al. Enhanced thermal decomposition of nitromethane on functionalized graphene sheets: ab initio molecular dynamics simulations[J].JournaloftheAmericanChemicalSociety, 2012, 134(46): 19011-19016.

[4] Zhang C Y, Wen Y S, Xue X G. Self-enhanced catalytic activities

of functionalized graphene sheets in the combustion of nitromethane: molecular dynamic simulations by molecular reactive force field[J].ACSAppliedMaterials&Interfaces, 2014, 6(15): 12235-12244.

[5] 劉英哲, 康瑩, 來蔚鵬, 等. 硝基甲烷在石墨烯表面初始反應機理的理論研究[J]. 含能材料, 2015, 23(9): 871-876.

LIU Ying-zhe, KANG Ying, LAI Wei-peng, et al. Reaction mechanism of nitromethane on the graphene surface: a theoretical study[J].ChineseJournalofEnergeticMaterials(HannengCailiao), 2015, 23(9): 871-876.

[6] Zhang C Y. Stress-induced activation of decomposition of organic explosives: a simple way to understand[J].JournalofMolecularModeling, 2013, 19(1): 477-483.

[7] Liu Y Z, Lai W P, Yu T, et al. Structural characteristics of liquid nitromethane at the nanoscale confinement in carbon nanotubes[J].JournalofMolecularModeling, 2014, 20(10): 2459.

[8] Vanommeslaeghe K, Hatcher E, Acharya C, et al. CHARMM general force field: a force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields[J].JournalofComputationalChemistry, 2010, 31(4): 671-690.

[9] Humphrey W, Dalke A, Schulten K. VMD: visual molecular dynamics[J].JournalofMolecularGraphics, 1996, 14(1): 33-38.

[10] MacKerell A D, Bashford D, Bellott M, et al. All-atom empirical potential for molecular modeling and dynamics studies of proteins[J].JournalofPhysicalChemistryB, 1998, 102(18): 3586-3616.

[11] Phillips J C, Braun R, Wang W, et al. Scalable molecular dynamics with NAMD[J].JournalofComputationalChemistry, 2005, 26(16): 1781-1802.

[12] Ryckaert J P, Ciccotti G, Berendsen H J C. Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics ofn-alkanes[J].JournalofComputationalPhysics, 1977, 23(3): 327-341.

[13] Andersen H C. Rattle: a “velocity” version of the shake algorithm for molecular dynamics calculations[J].JournalofComputationalPhysics, 1983, 52(1): 24-34.

[14] Darden T, York D, Pedersen L, et al. Particle mesh Ewald: anN·log(N) method for Ewald sums in large systems[J].JournalofChemicalPhysics, 1993, 98(12): 10089-10092.

[15] Tuckerman M, Berne B J, Martyna G J. Reversible multiple time scale molecular dynamics[J].JournalofChemicalPhysics, 1992, 97(3): 1990-2001.

[16] Feller S E, Zhang Y, Pastor R W, et al. Constant pressure molecular dynamics simulation: the langevin piston method[J].JournalofChemicalPhysics, 1995, 103(11): 4613-4621.

[17] Megyes T, Bálint S, Grósz T, et al. Structure of liquid nitromethane: comparison of simulation and diffraction studies[J].JournalofChemicalPhysics, 2007, 126(16): 164507.

[18] Fujiwara S, Sato T. Molecular dynamics simulations of structural formation of a single polymer chain: bond-orientational order and conformational defects[J].JournalofChemicalPhysics, 1997, 107(2): 613-622.