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

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

1 引言

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

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

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

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

2 計算模型和方法

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

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

TATB分子是一個苯的六個氫原子分別被硝基和氨基取代而形成的包含分子內(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)象的差異在于四個硝基的相對位置不同:α-HMX的四個硝基位于八元環(huán)的同一側(cè);β-HMX的四個硝基中相鄰的兩對分別位于八元環(huán)的兩側(cè);而γ-HMX的四個硝基中的一個位于八元環(huán)一側(cè),其余三個位于另一側(cè).

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

常溫常壓下TATB只有一種晶型,屬于三斜晶系,每個晶胞有2個分子,為P1空間群.17RDX分子在常溫常壓下形成一種穩(wěn)定的晶體,屬于正交晶系,每個晶胞有8個分子,為Pbca空間群.18而HMX分子能形成四種晶型(α-HMX、β-HMX、γ-HMX、δ-HMX),其中γ-HMX為水合晶型,19這里不做研究.HMX在常溫常壓下穩(wěn)定存在的是β型晶體,該晶體由β構(gòu)象的HMX分子構(gòu)成,每個晶胞2個分子,為單斜晶系P21/C空間群;20當(dāng)溫度升至375-377 K時,β-HMX晶體轉(zhuǎn)換成α-HMX晶體,它由α構(gòu)象構(gòu)成,每個晶胞含有8個α-HMX,為正交晶系Fdd2空間群;21當(dāng)溫度升至433-437 K時,α-HMX晶體轉(zhuǎn)換成δ-HMX晶體,它也由α構(gòu)象構(gòu)成,六方晶系P6122空間群,每個晶胞6個α-HMX分子.22γ構(gòu)象的HMX分子不形成晶體.分子模擬所研究的五種晶體結(jié)構(gòu)均引自文獻(xiàn),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 量子化學(xué)計算方法

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

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

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

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

計算范德華相互作用時,不同原子間的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計算的能量和能量導(dǎo)數(shù)確定分子內(nèi)的鍵參數(shù).電荷的初始參數(shù)用量子化學(xué)方法計算的Mulliken電荷,VDW的初始參數(shù)用OPLS力場29中相似的原子類型的非鍵參數(shù).這樣得到一個完整的力場后再用分子動力學(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 分子動力學(xué)模擬

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

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

3.1 原子類型

原子類型的定義僅僅與其附近的拓?fù)浣Y(jié)構(gòu)有關(guān),根據(jù)原子在分子中所處的環(huán)境定義.原子類型的定義由三部分組成:中心原子的元素符號、中心原子連接其它不同原子的個數(shù)和備注項.例如:c_4h2,c表示中心原子為碳,4表示碳原子一共連4個其它原子,h2表示碳原子所連的4個其它原子中有兩個是氫原子,所以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é)計算得到的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á)分子間的相互作用.因此電荷分布可以在一個較寬的范圍里采用,而分子間相互作用的其余部分由VDW勢來表達(dá).這一現(xiàn)象在OPLS29和COMPASS30的力場開發(fā)中已有報道.我們在優(yōu)化電荷和VDW參數(shù)時也看到由Mulliken電荷推導(dǎo)出的電荷參數(shù)過分高估分子間相互作用,以致無法調(diào)整VDW參數(shù)來得到合理的晶格能.這是因?yàn)樵谛纬煞肿泳w時分子間排列緊密,而這些分子都含強(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)計算.電荷用鍵增量(bond increment)30表達(dá).對每個原子,其所帶的電荷由式(3)計算得到,

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

3.3 分子性質(zhì)

力場中的分子內(nèi)參數(shù)是在分子間參數(shù)確定后通過擬合DFT計算的能量和能量導(dǎo)數(shù)確定的,完整的分子內(nèi)參數(shù)在Supporting Information中給出.在得到參數(shù)以后,用分子力學(xué)方法對TATB、RDX和HMX分子進(jìn)行結(jié)構(gòu)優(yōu)化并同量子化學(xué)計算結(jié)果對比來驗(yàn)證參數(shù)的可靠性.根據(jù)文獻(xiàn)報道,B3LYP/6-31G(d,p)計算可以在實(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ù)計算所得到的鍵長、鍵角、二面角和基態(tài)簡正振動頻率等結(jié)果與量子化學(xué)結(jié)果對比如圖2所示.力場優(yōu)化的鍵長與量子化學(xué)優(yōu)化的鍵長的均方差(MSD)小于0.0015 nm,鍵角的均方差小于3.7°,二面角均方差小于21.5°,振動頻率均方差小于95 cm-1.二面角較大的均方差源于在HMX的三種構(gòu)象中個別含H原子的二面角最大偏差達(dá)到40°左右.通過提高二面角的力常數(shù)可以將此偏差進(jìn)一步減小,但剛性太大的二面角函數(shù)影響了力場在描述凝聚相物理性質(zhì)時的表現(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晶體密度的計算結(jié)果為1.801 g·cm-3,比實(shí)驗(yàn)值偏低0.28%.β-HMX的密度比實(shí)驗(yàn)值低1.7%.總的來說,這些數(shù)據(jù)優(yōu)于文獻(xiàn)中報道的力場8,11,33得到的結(jié)果,只有Gee力場12計算得到了和實(shí)驗(yàn)完全一致的數(shù)據(jù),是例外.對于模擬溫度高于常溫的兩種材料α-HMX和δ-HMX,密度相比于實(shí)驗(yàn)值較大,分別是偏小5.6%和3.6%.

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

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

上式中,Einter是體系的內(nèi)聚能,即分子間相互作用;R為摩爾氣體常數(shù),8.314 J·K-1·mol-1;T為溫度.α-HMX和δ-HMX的升華焓缺乏準(zhǔn)確的實(shí)驗(yàn)測量值.TATB、RDX和β-HMX的升華焓計算結(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ì)

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

含能材料晶體在高壓下的熱力學(xué)狀態(tài)方程是描述該類材料的一個重要性質(zhì).我們計算了298 K下幾種含能材料分子晶體的壓力-體積(p-V)曲線.圖3是計算得到的TATB、RDX和β-HMX的p-V曲線及其和文獻(xiàn)35-39的對比.其中RDX在大于4 GPa的壓力下會發(fā)生相變,因此只計算了小于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)中報道的計算值.

圖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)的理論計算值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曲線所得的體積模量的計算結(jié)果與文獻(xiàn)報道結(jié)果35-39的對比.由表中數(shù)據(jù)可以看出,該力場在預(yù)測材料機(jī)械性質(zhì)上有著較好的表現(xiàn),所得結(jié)果和實(shí)驗(yàn)值以及其它理論計算值相比都比較接近.

5 結(jié)論

報道了一個適用于TATB、RDX和HMX含能材料的,可以在通用軟件中使用的分子力學(xué)力場,并通過計算分子和分子晶體的物理性質(zhì)驗(yàn)證了該力場.驗(yàn)證結(jié)果顯示該力場可以較好地重現(xiàn)量子化學(xué)DFT預(yù)測的分子結(jié)構(gòu)、構(gòu)象和振動頻率,得到和實(shí)驗(yàn)值基本一致的(包括HMX在不同溫度下的三種晶型)的晶胞參數(shù)、密度和升華焓.擴(kuò)展到預(yù)測TATB、RDX和β-HMX分子晶體的等溫p-V曲線、體積模量及其對壓力的一階導(dǎo)也得到與文獻(xiàn)報道的計算和實(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é)報,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é)報,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