適用于TATB,RDX,HMX 含能材料的全原子力場的建立與驗證

金 釗 劉 建 王麗莉 曹風雷 孫 淮,*

(1上海交通大學化學化工學院,上海200240;2中國工程物理研究院計算機應用研究所,四川綿陽621900)

1 引言

含能材料在航天、國防、能源領域有著重要的應用.采用理論與計算化學方法研究含能材料能夠彌補實驗研究方法的不足,提高研發效率.與量子化學方法對比,基于分子力學力場的模擬方法能夠在較大尺度下預測材料特別是含能材料與高分子復合材料的物理性質,1-7但力場參數的不足嚴重阻礙了這一類方法的應用.

文獻中已有一些關于含能材料的分子力學力場的報道.Sorescu等8,9采用剛性分子模型開發了環三亞甲基三硝胺(RDX)和環四亞甲基四硝胺(HMX)的力場,分子間的相互作用使用庫侖和Buckingham勢函數(簡稱SRT力場)來表示.計算結果顯示該力場可以比較準確地描述晶體的結構.Smith和Bharadwaj10開發了HMX的全原子力場,其力場函數包括了分子內和分子間兩部分,其中分子內部分用簡諧振動函數描述鍵長、鍵角,余弦函數描述二面角,而分子間的非鍵項采用Buckingham勢函數和修正的庫侖項(簡稱SB力場).Boyd等11開發了RDX的全原子力場,分子內鍵參數采用Morse函數,角參數為簡諧振動函數,二面角采用余弦函數形式,分子間非鍵作用函數形式為Buckingham勢函數和庫侖項(簡稱Boyd力場).Gee等12開發了三硝基三氨基苯(TATB)的力場,其分子內勢函數形式和SB力場一樣,但分子間勢函數采用的是Lennard-Jones(12-6)勢函數和庫侖項(簡稱Gee力場),該力場在描述TATB晶體結構方面具有相當的穩定性和準確性.這些工作表明:用經驗的力場函數模型可以較好地描述這一類含能材料的物理性質.但是已報道的力場均是針對單一分子開發的,而使用的函數形式也不一致,因此這些力場不具備可遷移性和在通用軟件中的普適性.

為了擴大分子模擬方法在含能材料領域里的應用,需要建立具有可遷移性和普適性的分子力學力場.只有具有可遷移性,力場才可以用來預測不同條件下的物理性質.構建可遷移性力場的要點是同時擬合一組相似的分子,在不同的分子中根據相同的化學環境定義一致的原子類型,以獲得在不同的條件下均可適用的參數.本文以TATB、RDX、HMX為對象,開發一個對此類材料分子適用的分子力學力場.為了普適性的要求,我們采用常用的力場函數形式,這樣得到的力場可以在通用的分子模擬軟件(如LAMMPS,13,14NAMD15和GROMACS16)上直接使用.

本文以幾種常見的含能材料為研究對象,采用量子化學計算和分子動力學模擬等方法,開發了可遷移的力場.應用所得力場計算了分子和分子晶體的性質,并通過計算材料的p-V狀態方程和機械性質對該力場進一步驗證.

2 計算模型和方法

2.1 分子結構

圖1 TATB、RDX和HMX的分子結構Fig.1 Molecular structures of TATB,RDX,and HMX

TATB分子是一個苯的六個氫原子分別被硝基和氨基取代而形成的包含分子內氫鍵的平面結構,RDX和HMX是含有硝胺基團的六元環和八元環.分子的結構見圖1.TATB和RDX的穩定結構單一,但HMX有三種穩定的構象,分別為α-HMX、β-HMX和γ-HMX.這三種構象的差異在于四個硝基的相對位置不同:α-HMX的四個硝基位于八元環的同一側;β-HMX的四個硝基中相鄰的兩對分別位于八元環的兩側;而γ-HMX的四個硝基中的一個位于八元環一側,其余三個位于另一側.

2.2 晶體結構

常溫常壓下TATB只有一種晶型,屬于三斜晶系,每個晶胞有2個分子,為P1空間群.17RDX分子在常溫常壓下形成一種穩定的晶體,屬于正交晶系,每個晶胞有8個分子,為Pbca空間群.18而HMX分子能形成四種晶型(α-HMX、β-HMX、γ-HMX、δ-HMX),其中γ-HMX為水合晶型,19這里不做研究.HMX在常溫常壓下穩定存在的是β型晶體,該晶體由β構象的HMX分子構成,每個晶胞2個分子,為單斜晶系P21/C空間群;20當溫度升至375-377 K時,β-HMX晶體轉換成α-HMX晶體,它由α構象構成,每個晶胞含有8個α-HMX,為正交晶系Fdd2空間群;21當溫度升至433-437 K時,α-HMX晶體轉換成δ-HMX晶體,它也由α構象構成,六方晶系P6122空間群,每個晶胞6個α-HMX分子.22γ構象的HMX分子不形成晶體.分子模擬所研究的五種晶體結構均引自文獻,17,18,20-22建模使用Direct Force Field 7.1軟件包.23超晶胞的組成如下:TATB含有3×3×4個單胞、RDX 含有2×2×2個單胞、β-HMX含有4×2×3個單胞、α-HMX含有2×1×4個單胞、δ-HMX含有3×3×1個單胞.

2.3 量子化學計算方法

量子化學計算采用B3LYP密度泛函方法和6-31G(d,p)基組24,25優化單分子結構、計算ESP26和Mulliken27電荷分布、構象能及分子的總能量對笛卡爾坐標的一階和二階導數.這些數據用來擬合力場的鍵參數和電荷參數.量子化學計算采用Gaussian 03軟件包28實施.

2.4 力場勢函數形式

采用的勢函數包括成鍵項(分子內)和非鍵項(分子間)兩部分.成鍵項包括鍵、角、二面角、面外鍵角和交叉項等.非鍵項包括靜電相互作用(庫侖形式)和范德華相互作用(LJ 12-6形式):

式中鍵參數部分kb、ka、kt、ko、kbb、kba分別為鍵伸縮、角伸縮、二面角扭轉、面外鍵角、鍵鍵交叉項和鍵角交叉項的力常數,b為鍵長,θ為鍵角,φ為二面角,χ為面外鍵角,下標0表示相應的平衡值,c、d、e、f為系數;非鍵參數部分q為原子所帶電荷數;rij是非鍵原子間的距離,r0ij是原子i和原子j的范德華相互作用半徑,εij是原子i和原子j的相互作用勢阱.

計算范德華相互作用時,不同原子間的Lennard-Jones參數采用Lorentz-Berthlot混合規則得到:

2.5 參數化方法

力場參數化是對分子內和分子間的參數分別優化并迭代的過程.為了保持力場參數的一致性,亞甲基和氨基的參數由我們以前推導的烷烴和胺類化合物的參數遷移得到.其它參數通過如下步驟得到:首先,確定并固定初始的電荷及范德華(VDW)參數,用最小二乘法擬合量子化學DFT計算的能量和能量導數確定分子內的鍵參數.電荷的初始參數用量子化學方法計算的Mulliken電荷,VDW的初始參數用OPLS力場29中相似的原子類型的非鍵參數.這樣得到一個完整的力場后再用分子動力學模擬TATB、RDX和HMX晶體的晶胞性質、升華焓,并同實驗數據對比優化電荷和VDW參數.上述過程循環迭代得到優化的力場參數.擬合過程使用Direct Force Field 7.1軟件包.

2.6 分子動力學模擬

晶體的分子動力學模擬部分采用LAMMPS軟件包,長程電荷作用采用Ewald/n方法,截斷半徑(cut off)為1.0 nm,時間步長1 fs,模擬時間2 ns,其中1 ns平衡,1 ns采樣.采用恒溫恒壓系綜(NPT),用Nose/Hoover控溫控壓方法,允許超晶胞的6個自由度(即3個邊長和3個夾角)變化.

3 力場的建立和初步驗證

3.1 原子類型

原子類型的定義僅僅與其附近的拓撲結構有關,根據原子在分子中所處的環境定義.原子類型的定義由三部分組成:中心原子的元素符號、中心原子連接其它不同原子的個數和備注項.例如:c_4h2,c表示中心原子為碳,4表示碳原子一共連4個其它原子,h2表示碳原子所連的4個其它原子中有兩個是氫原子,所以c_4h2表示的是亞甲基上的碳原子類型.本文涉及到的原子類型有如下幾種:苯環上連氨基的碳原子c_3an,苯環上連硝基的碳原子c_3ani,亞甲基上碳原子c_4h2,亞甲基氫原子h_1,氨基氫原子h_1n,氨基氮原子n_3h2,連硝基的氮原子n_3no,硝基上氮原子n_3o,硝基上氧原子o_1n.

3.2 電荷和VDW參數的確定

量子化學計算得到的ESP和Mulliken電荷見表1.從數據可以看出,對于相同的原子類型,Mulliken電荷比ESP電荷在不同的分子間更穩定.這是由于ESP電荷是通過擬合靜電勢得到的,擬合的結果受到采樣方式的影響.26而根據波函數確定的Mulliken電荷更能準確地反映原子間電子密度的分布,因此我們選擇Mulliken電荷作為初始電荷.

分子間相互作用本質上都是源于庫侖力,可以用電荷-電荷、電荷-偶極、偶極-偶極以及包括四極矩和誘導偶極在內的庫侖力來表達.而在力場方法中僅用到了位于原子上的電荷和VDW勢、電荷分布和VDW勢的結合用來表達分子間的相互作用.因此電荷分布可以在一個較寬的范圍里采用,而分子間相互作用的其余部分由VDW勢來表達.這一現象在OPLS29和COMPASS30的力場開發中已有報道.我們在優化電荷和VDW參數時也看到由Mulliken電荷推導出的電荷參數過分高估分子間相互作用,以致無法調整VDW參數來得到合理的晶格能.這是因為在形成分子晶體時分子間排列緊密,而這些分子都含強極性基團硝基,點電荷模型在近距離產生過強的靜電作用.我們把點電荷參數統一地縮小至0.65倍,再調整VDW參數,得到了較好的效果.

表2列出優化的VDW勢(即LJ 12-6函數)的參數和電荷參數.為了參數的完整性,從烷烴和胺類遷移的參數也一并列出.LJ 12-6參數按原子類型給出,不同原子類型間的相互作用采用組合規則(式2)計算.電荷用鍵增量(bond increment)30表達.對每個原子,其所帶的電荷由式(3)計算得到,

式中j表示與i相連的所有其它原子.

3.3 分子性質

力場中的分子內參數是在分子間參數確定后通過擬合DFT計算的能量和能量導數確定的,完整的分子內參數在Supporting Information中給出.在得到參數以后,用分子力學方法對TATB、RDX和HMX分子進行結構優化并同量子化學計算結果對比來驗證參數的可靠性.根據文獻報道,B3LYP/6-31G(d,p)計算可以在實驗精度范圍里較好地再現分子的結構和構象能的實驗值.24,25

表1 給定原子類型在不同分子里的ESP和Mulliken電荷Table 1 ESPand Mulliken charges for atom types in different molecules

表2 范德華參數和電荷鍵增量參數Table 2 VDW and bond increment(BI)parameters of all the atom types

用優化參數計算所得到的鍵長、鍵角、二面角和基態簡正振動頻率等結果與量子化學結果對比如圖2所示.力場優化的鍵長與量子化學優化的鍵長的均方差(MSD)小于0.0015 nm,鍵角的均方差小于3.7°,二面角均方差小于21.5°,振動頻率均方差小于95 cm-1.二面角較大的均方差源于在HMX的三種構象中個別含H原子的二面角最大偏差達到40°左右.通過提高二面角的力常數可以將此偏差進一步減小,但剛性太大的二面角函數影響了力場在描述凝聚相物理性質時的表現.這一問題也反映了力場方法的局限性.

3.4 晶胞結構和升華焓

晶胞結構和升華焓是用來優化非鍵參數的指標.用優化的參數在實驗溫度21,22和常壓下對TATB、RDX和HMX晶體進行模擬得到的晶胞參數、晶體密度及升華焓在表3列出,并與實驗值17,18,20-22,31,32進行對比.

由結果可以看出,力場能準確描述這些晶體的結構性質,模擬的晶胞邊長結果與實驗值差別全部小于3.6%.晶胞夾角最大差異為TATB的β角,與實驗測量值相差4.1°,其余夾角偏差全部小于2.3°.在室溫下得到的晶體密度與實驗值相吻合.TATB的密度為1.923 g·cm-3,比實驗值偏低0.72%.RDX晶體密度的計算結果為1.801 g·cm-3,比實驗值偏低0.28%.β-HMX的密度比實驗值低1.7%.總的來說,這些數據優于文獻中報道的力場8,11,33得到的結果,只有Gee力場12計算得到了和實驗完全一致的數據,是例外.對于模擬溫度高于常溫的兩種材料α-HMX和δ-HMX,密度相比于實驗值較大,分別是偏小5.6%和3.6%.

圖2 用力場(FF)和量子化學(QM)方法優化得到的分子鍵長(a),鍵角(b),二面角(c)及振動頻率(d)的對比Fig.2 Comparisons of optimized bond lengths(a),bond angles(b),dihedral angles(c),and vibrational frequencies(d)between quantum chemistry(QM)and force field(FF)methods

表3 晶體的晶胞參數、密度和升華焓(ΔHsub)Table 3 Crystal cell parameters,densities,heat of sublimation(ΔHsub)

表3中升華焓ΔHsub按照如下公式進行近似計算:

上式中,Einter是體系的內聚能,即分子間相互作用;R為摩爾氣體常數,8.314 J·K-1·mol-1;T為溫度.α-HMX和δ-HMX的升華焓缺乏準確的實驗測量值.TATB、RDX和β-HMX的升華焓計算結果和實驗值相比偏差分別為2.9%、1.2%和-4.1%.而之前的力場預測的誤差分別是TATB:-22.7%(Gee34),RDX:-11.4%(SRT8),-10.9%(Boyd11)和β-HMX:3.0%(SB33).

4 晶體的狀態方程和機械性質

我們通過計算晶體的狀態方程和機械性質進一步驗證了得到的力場的可靠性.這些性質的計算對優化參數有指導意義,但沒有直接用來優化參數,因此提供了更為嚴格的測試條件.

含能材料晶體在高壓下的熱力學狀態方程是描述該類材料的一個重要性質.我們計算了298 K下幾種含能材料分子晶體的壓力-體積(p-V)曲線.圖3是計算得到的TATB、RDX和β-HMX的p-V曲線及其和文獻35-39的對比.其中RDX在大于4 GPa的壓力下會發生相變,因此只計算了小于4 GPa的情況.從圖中可以看出,TATB的預測結果和實驗值基本一致;而RDX在高壓下和實驗值稍有偏差,高估體積約0.9%;β-HMX的結果和實驗值相比在低壓下(0.5 GPa)體積偏小0.8%左右.總的來說,這些數據優于文獻中報道的計算值.

圖3 TATB、RDX和β-HMX在298 K下的p-V曲線Fig.3 p-V curves of TATB,RDX,and β-HMX at 298 K

表4 298 K下晶體的體積模量(B0)及其一階導數(B′0)的理論計算值Table 4 Calculated bulk modulus(B0)and their first derivatives(B′0)at 298 K

體積模量及其對壓力的一階導數利用三階Birch-Murnaghan公式擬合p-V曲線得到:式中B0為晶體的體積模量,B'0為體積模量對壓力的一階導數,p為壓力,V為該壓力下材料的體積,V0為初始條件下材料的體積.表4為擬合三種材料的p-V曲線所得的體積模量的計算結果與文獻報道結果35-39的對比.由表中數據可以看出,該力場在預測材料機械性質上有著較好的表現,所得結果和實驗值以及其它理論計算值相比都比較接近.

5 結論

報道了一個適用于TATB、RDX和HMX含能材料的,可以在通用軟件中使用的分子力學力場,并通過計算分子和分子晶體的物理性質驗證了該力場.驗證結果顯示該力場可以較好地重現量子化學DFT預測的分子結構、構象和振動頻率,得到和實驗值基本一致的(包括HMX在不同溫度下的三種晶型)的晶胞參數、密度和升華焓.擴展到預測TATB、RDX和β-HMX分子晶體的等溫p-V曲線、體積模量及其對壓力的一階導也得到與文獻報道的計算和實驗數據一致的結果.這些結果顯示了用經驗的力場函數可以在較大溫度壓力范圍內描述含能材料的物理性質.和文獻中已發表的力場對比,由于采用了統一的力場函數形式和原子類型描述含硝基、氨基的芳香或非芳香性環狀化合物,本文提出的力場具有更好的可遷移性和普適性,為我們進一步深入研究這些含能材料的物理性質提供了可能,也為我們進一步開發粗粒化力場、在更大尺度上研究這一類材料打下了基礎.

Supporting Information: The valence parameters of the force field have been included.This information is available free of charge via the internet at http://www.whxb.pku.edu.cn.

(1)Gee,R.H.;Maiti,A.;Bastea,S.;Fried,L.E.Macromolecules 2007,40,3422.doi:10.1021/ma0702501

(2) Maiti,A.;Gee,R.H.;Hoffman,D.M.;Fried,L.E.J.Appl.Phys.2008,103,053504.doi:10.1063/1.2838319

(3) Xiao,J.;Huang,H.;Li,J.;Zhang,H.;Zhu,W.;Xiao,H.J.Mater.Sci.2008,43,5685.doi:10.1007/s10853-008-2704-0

(4) Xu,X.;Xiao,J.;Huang,H.;Li,J.;Xiao,H.J.Hazard.Mater.2010,175,423.doi:10.1016/j.jhazmat.2009.10.023

(5) Zhou,Y.;Long,X.;Wei,X.J.Mol.Model.2011,17,3015.doi:10.1007/s00894-011-0977-8

(6) Huang,Y.C.;Hu,Y.J.;Xiao,J.J.;Yin,K.L.;Xiao,H.M.Acta Phys.-Chim.Sin.2005,21,425.[黃玉成,胡應杰,肖繼軍,殷開梁,肖鶴鳴.物理化學學報,2005,21,425.]doi:10.3866/PKU.WHXB20050416

(7) Xiao,J.J.;Gu,C.G.;Fang,G.Y.;Zhu,W.;Xiao,H.M.Acta Chim.Sin.2005,63,439.[肖繼軍,谷成剛,方國勇,朱 偉,肖鶴鳴.化學學報,2005,63,439.]

(8) Sorescu,D.C.;Rice,B.M.;Thompson,D.L.J.Phys.Chem.B 1997,101,798.doi:10.1021/jp9624865

(9) Sorescu,D.C.;Rice,B.M.;Thompson,D.L.J.Phys.Chem.B 1998,102,6692.doi:10.1021/jp981661+

(10)Smith,G.D.;Bharadwaj,R.K.J.Phys.Chem.B 1999,103,3570.doi:10.1021/jp984599p

(11) Boyd,S.;Gravelle,M.;Politzer,P.J.Chem.Phys.2006,124,104508.doi:10.1063/1.2176621

(12) Gee,R.H.;Roszak,S.;Balasubramanian,K.;Fried,L.E.J.Chem.Phys.2004,120,7059.doi:10.1063/1.1676120

(13) Plimpton,S.J.Comput.Phys.1995,117,1.doi:10.1006/jcph.1995.1039

(14)Plimpton,S.LAMMPS:Large-ScaleAtomic/Molecular Massively Parallel Simulator.http://lammps.sandia.gov.

(15) Phillips,J.C.;Braun,R.;Wang,W.;Gumbart,J.;Tajkhorshid,E.;Villa,E.;Chipot,C.;Skeel,R.D.;Kale,L.;Schulten,K.J.Comput.Chem.2005,26,1781.

(16) Van Der Spoel,D.;Lindahl,E.;Hess,B.;Groenhof,G.;Mark,A.E.;Berendsen,H.J.J.Comput.Chem.2005,26,1701.

(17) Cady,H.H.;Larson,A.C.Acta Crystallogr.1965,18,485.doi:10.1107/S0365110X6500107X

(18) Choi,C.;Prince,E.Acta Crystallogr.Sect.B 1972,28,2857.doi:10.1107/S0567740872007046

(19) Cady,H.H.;Smith,L.C.Los Alamos Scientific Laboratory Report LAMS-2652 TID-4500;LosAlamos National Laboratory:LosAlamos,NM,1961.

(20) Choi,C.S.;Boutin,H.P.Acta Crystallogr.Sect.B 1970,26,1235.doi:10.1107/S0567740870003941

(21) Cady,H.H.;Larson,A.C.;Cromer,D.T.Acta Crystallogr.1963,16,617.doi:10.1107/S0365110X63001651

(22) Cobbledick,R.;Small,R.Acta Crystallogr.Sect.B 1974,30,1918.doi:10.1107/S056774087400611X

(23) Direct Force Field,7.1;Aeon Technology Inc.:San Diego,CA,USA,2012.

(24) Rice,B.M.;Chabalowski,C.F.J.Phys.Chem.A 1997,101,8720.doi:10.1021/jp972062q

(25)Riley,K.E.;Op't Holt,B.T.;Merz,K.M.J.Chem.Theory Comput.2007,3,407.doi:10.1021/ct600185a

(26)Breneman,C.M.;Wiberg,K.B.J.Comput.Chem.1990,11,361.

(27) Mulliken,R.S.J.Chem.Phys.1955,23,1833.doi:10.1063/1.1740588

(28) Frisch,M.;Trucks,G.;Schlegel,H.;et al.Gaussian 03,Revision C.02;Gaussian Inc.:Wallingford,CT,2004.

(29)Jorgensen,W.L.;Maxwell,D.S.;Tirado-Rives,J.J.Am.Chem.Soc.1996,118,11225.doi:10.1021/ja9621760

(30) Sun,H.J.Phys.Chem.B 1998,102,7338.doi:10.1021/jp980939v

(31) Rosen,J.M.;Dickinson,C.J.Chem.Eng.Data 1969,14,120.doi:10.1021/je60040a044

(32) Taylor,J.W.;Crookes,R.J.J.Chem.Soc.Faraday Trans.1976,72,723.doi:10.1039/f19767200723

(33) Bedrov,D.;Ayyagari,C.;Smith,G.D.;Sewell,T.D.;Menikoff,R.;Zaug,J.M.J.Comput.Aided Mater.Des.2001,8,77.doi:10.1023/A:1020046817543

(34) Rai,N.;Bhatt,D.;Siepmann,J.I.;Fried,L.E.J.Chem.Phys.2008,129,194510.doi:10.1063/1.3006054

(35) Stevens,L.L.;Velisavljevic,N.;Hooks,D.E.;Dattelbaum,D.M.Propellants Explos.Pyrotech.2008,33,286.doi:10.1002/prep.v33:4

(36) Bedrov,D.;Borodin,O.;Smith,G.D.;Sewell,T.D.;Dattelbaum,D.M.;Stevens,L.L.J.Chem.Phys.2009,131,224703.doi:10.1063/1.3264972

(37) Olinger,B.;Roof,B.;Cady,H.In Proceedings of International Symposium on High Dynamic Pressures;Commissariat a l′EnergieAtomique:Paris,France,1978;p 3.

(38) Zheng,L.;Thompson,D.L.J.Chem.Phys.2006,125,084505.doi:10.1063/1.2238860

(39)Sorescu,D.C.;Rice,B.M.;Thompson,D.L.J.Phys.Chem.B 1999,103,6783.doi:10.1021/jp991202o