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

金 釗 劉 建 王麗莉 曹風(fēng)雷 孫 淮,*

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

1 引言

含能材料在航天、國防、能源領(lǐng)域有著重要的應(yīng)用.采用理論與計(jì)算化學(xué)方法研究含能材料能夠彌補(bǔ)實(shí)驗(yàn)研究方法的不足,提高研發(fā)效率.與量子化學(xué)方法對比,基于分子力學(xué)力場的模擬方法能夠在較大尺度下預(yù)測材料特別是含能材料與高分子復(fù)合材料的物理性質(zhì),1-7但力場參數(shù)的不足嚴(yán)重阻礙了這一類方法的應(yīng)用.

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

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

本文以幾種常見的含能材料為研究對象,采用量子化學(xué)計(jì)算和分子動(dòng)力學(xué)模擬等方法,開發(fā)了可遷移的力場.應(yīng)用所得力場計(jì)算了分子和分子晶體的性質(zhì),并通過計(jì)算材料的p-V狀態(tài)方程和機(jī)械性質(zhì)對該力場進(jìn)一步驗(yàn)證.

2 計(jì)算模型和方法

2.1 分子結(jié)構(gòu)

圖1 TATB、RDX和HMX的分子結(jié)構(gòu)Fig.1 Molecular structures of TATB,RDX,and HMX

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

2.2 晶體結(jié)構(gòu)

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

2.3 量子化學(xué)計(jì)算方法

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

2.4 力場勢函數(shù)形式

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

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

計(jì)算范德華相互作用時(shí),不同原子間的Lennard-Jones參數(shù)采用Lorentz-Berthlot混合規(guī)則得到:

2.5 參數(shù)化方法

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

2.6 分子動(dòng)力學(xué)模擬

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

3 力場的建立和初步驗(yàn)證

3.1 原子類型

原子類型的定義僅僅與其附近的拓?fù)浣Y(jié)構(gòu)有關(guān),根據(jù)原子在分子中所處的環(huán)境定義.原子類型的定義由三部分組成:中心原子的元素符號、中心原子連接其它不同原子的個(gè)數(shù)和備注項(xiàng).例如:c_4h2,c表示中心原子為碳,4表示碳原子一共連4個(gè)其它原子,h2表示碳原子所連的4個(gè)其它原子中有兩個(gè)是氫原子,所以c_4h2表示的是亞甲基上的碳原子類型.本文涉及到的原子類型有如下幾種:苯環(huán)上連氨基的碳原子c_3an,苯環(huán)上連硝基的碳原子c_3ani,亞甲基上碳原子c_4h2,亞甲基氫原子h_1,氨基氫原子h_1n,氨基氮原子n_3h2,連硝基的氮原子n_3no,硝基上氮原子n_3o,硝基上氧原子o_1n.

3.2 電荷和VDW參數(shù)的確定

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

分子間相互作用本質(zhì)上都是源于庫侖力,可以用電荷-電荷、電荷-偶極、偶極-偶極以及包括四極矩和誘導(dǎo)偶極在內(nèi)的庫侖力來表達(dá).而在力場方法中僅用到了位于原子上的電荷和VDW勢、電荷分布和VDW勢的結(jié)合用來表達(dá)分子間的相互作用.因此電荷分布可以在一個(gè)較寬的范圍里采用,而分子間相互作用的其余部分由VDW勢來表達(dá).這一現(xiàn)象在OPLS29和COMPASS30的力場開發(fā)中已有報(bào)道.我們在優(yōu)化電荷和VDW參數(shù)時(shí)也看到由Mulliken電荷推導(dǎo)出的電荷參數(shù)過分高估分子間相互作用,以致無法調(diào)整VDW參數(shù)來得到合理的晶格能.這是因?yàn)樵谛纬煞肿泳w時(shí)分子間排列緊密,而這些分子都含強(qiáng)極性基團(tuán)硝基,點(diǎn)電荷模型在近距離產(chǎn)生過強(qiáng)的靜電作用.我們把點(diǎn)電荷參數(shù)統(tǒng)一地縮小至0.65倍,再調(diào)整VDW參數(shù),得到了較好的效果.

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

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

3.3 分子性質(zhì)

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

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

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

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

3.4 晶胞結(jié)構(gòu)和升華焓

晶胞結(jié)構(gòu)和升華焓是用來優(yōu)化非鍵參數(shù)的指標(biāo).用優(yōu)化的參數(shù)在實(shí)驗(yàn)溫度21,22和常壓下對TATB、RDX和HMX晶體進(jìn)行模擬得到的晶胞參數(shù)、晶體密度及升華焓在表3列出,并與實(shí)驗(yàn)值17,18,20-22,31,32進(jìn)行對比.

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

圖2 用力場(FF)和量子化學(xué)(QM)方法優(yōu)化得到的分子鍵長(a),鍵角(b),二面角(c)及振動(dòng)頻率(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 晶體的晶胞參數(shù)、密度和升華焓(ΔHsub)Table 3 Crystal cell parameters,densities,heat of sublimation(ΔHsub)

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

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

4 晶體的狀態(tài)方程和機(jī)械性質(zhì)

我們通過計(jì)算晶體的狀態(tài)方程和機(jī)械性質(zhì)進(jìn)一步驗(yàn)證了得到的力場的可靠性.這些性質(zhì)的計(jì)算對優(yōu)化參數(shù)有指導(dǎo)意義,但沒有直接用來優(yōu)化參數(shù),因此提供了更為嚴(yán)格的測試條件.

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

圖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)及其一階導(dǎo)數(shù)(B′0)的理論計(jì)算值Table 4 Calculated bulk modulus(B0)and their first derivatives(B′0)at 298 K

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

5 結(jié)論

報(bào)道了一個(gè)適用于TATB、RDX和HMX含能材料的,可以在通用軟件中使用的分子力學(xué)力場,并通過計(jì)算分子和分子晶體的物理性質(zhì)驗(yàn)證了該力場.驗(yàn)證結(jié)果顯示該力場可以較好地重現(xiàn)量子化學(xué)DFT預(yù)測的分子結(jié)構(gòu)、構(gòu)象和振動(dòng)頻率,得到和實(shí)驗(yàn)值基本一致的(包括HMX在不同溫度下的三種晶型)的晶胞參數(shù)、密度和升華焓.擴(kuò)展到預(yù)測TATB、RDX和β-HMX分子晶體的等溫p-V曲線、體積模量及其對壓力的一階導(dǎo)也得到與文獻(xiàn)報(bào)道的計(jì)算和實(shí)驗(yàn)數(shù)據(jù)一致的結(jié)果.這些結(jié)果顯示了用經(jīng)驗(yàn)的力場函數(shù)可以在較大溫度壓力范圍內(nèi)描述含能材料的物理性質(zhì).和文獻(xiàn)中已發(fā)表的力場對比,由于采用了統(tǒng)一的力場函數(shù)形式和原子類型描述含硝基、氨基的芳香或非芳香性環(huán)狀化合物,本文提出的力場具有更好的可遷移性和普適性,為我們進(jìn)一步深入研究這些含能材料的物理性質(zhì)提供了可能,也為我們進(jìn)一步開發(fā)粗?；觥⒃诟蟪叨壬涎芯窟@一類材料打下了基礎(chǔ).

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.[黃玉成,胡應(yīng)杰,肖繼軍,殷開梁,肖鶴鳴.物理化學(xué)學(xué)報(bào),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.[肖繼軍,谷成剛,方國勇,朱 偉,肖鶴鳴.化學(xué)學(xué)報(bào),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