材料的計算機模擬方法課件

2022/11/141材料的計算機模擬方法姬廣富沖擊波物理與爆轟物理國防科技重點實驗室內(nèi)容材料計算方法以及最新的進展；含能材料的引發(fā)反應(yīng)；靜高壓下材料的行為；材料的降解動力學(xué)模擬；氣態(tài)CO在高溫高壓下的歧化反應(yīng)；

HMX固體的性能理論模擬1計算方法簡介在含能材料的模擬計算中，目前主要有三種方法：第一性原理（或從頭算）的電子結(jié)構(gòu)計算(abinitio)、半經(jīng)驗方法分子力學(xué)方法(MM)。

通過對材料的模擬計算可以得到如下的性質(zhì)：生成熱、鍵能和反應(yīng)能、分子（或晶體）的能量和結(jié)構(gòu)熱化學(xué)穩(wěn)定性)、過渡態(tài)的能量和結(jié)構(gòu)（活化能）、反應(yīng)路徑以及動力學(xué)、機理、分子中的電荷分布（反應(yīng)點）、取代效應(yīng)、振動頻率（紅外和Raman光譜）、電子的激發(fā)躍遷（UV/可見光譜）、磁屏蔽效應(yīng)（NMR譜）。。。。。。通過對這些性質(zhì)的預(yù)測，在含能材料中有許多實用價值，如：可以研究合成路徑、反應(yīng)產(chǎn)物、爆轟產(chǎn)物狀態(tài)方程、爆轟引發(fā)的機理等等。1.1從頭算(abinitio)方法從頭算方法在計算中是基于量子力學(xué)基本定律，沒有引入經(jīng)驗或?qū)嶒瀰?shù)，僅使用5個常數(shù)（電子質(zhì)量：9.10938215E-31Kg、電子電量：1.602176487E-19C、光速：c=2.99792458×108m/s、Planck常數(shù)：h=6.62606896×10-34J·s和波爾茲曼常數(shù)1.3806504E-23J/K），對化學(xué)體系利用滿足嚴(yán)格數(shù)學(xué)近似的基組來解Schr?dinger方程：動能是所有粒子的動能加和勢能代表粒子之間的庫侖相互作用（即：核與核之間、電子與電子之間的排斥以及電子與核之間的吸引）對電子而言，q=-e，對核而言，原子序數(shù)為Z：q=+Ze。通過對Schr?dinger方程的形式觀察，我們知道對于多于兩個以上粒子的分子，嚴(yán)格解Schr?dinger方程是不可能的，必須引入一些必要的合理近似，甚至一些經(jīng)驗參數(shù)，如非相對論解、Born-Oppenheimer近似、單粒子近似、原子軌道的線性組合等。對解出的Schr?dinger方程是否合適，必須滿足以下判據(jù)：這個解必須很好地被定義且能夠用來說明分子結(jié)構(gòu)和分子中電子的狀態(tài)。分子的勢能按照原子核的位移必須平滑連續(xù)可變。模型不能有偏離（即：假設(shè)化學(xué)鍵只存在于原子之間）。模型必須與體系大小相關(guān)（即解和相關(guān)誤差必須與分子的大小成比例），模型必須是可變的（即：近似解必須提供體系能量真值的上界，因此近似解應(yīng)具有最低能量值與真實的波函數(shù)相一致）。1.1.1基組一般來說，基組屬于數(shù)學(xué)函數(shù)的一類，用來解微分方程?；M：G、PW、G+PW量化計算，基組一般約定為高斯基組來代表原子軌道，用來優(yōu)化和產(chǎn)生所需體系的化學(xué)性質(zhì)物理計算（DFT方法）,特別是金屬的計算，經(jīng)常用PW；GPW標(biāo)準(zhǔn)的從頭算軟件包都提供基組、當(dāng)分子中包含較大的原子（Z>30）時，我們經(jīng)常將原子的內(nèi)層電子合并作近似處理（贗勢）。帶極化函數(shù)的基組，在基組中加入極化函數(shù)后可以允許原子核外電荷的不均勻位移，因此它可以改進對化學(xué)鍵的描述。極化函數(shù)可以描述那些和孤立原子相比更高的角動量量子數(shù)的軌道（如p-型函數(shù)用于H和He，d-型函數(shù)用于Z>2的原子），并將之加入到價電子殼層中。例如，6-31G(d)基組的構(gòu)建是將6個d-型簡單高斯函數(shù)加入到6-31G中來描述每個非氫原子；同樣地，6-31G(d,p)是6-31G(d)對重原子,但加入了一組p-型高斯函數(shù)到氫原子或氦原子。將p-軌道加入到氫原子中對那些含有以氫原子為橋原子的體系而言特別有用。帶彌散函數(shù)的基組，對于那些電子密度明顯遠離核的中心的原子或離子（如陰離子、孤電子對以及激發(fā)態(tài)），就需要用彌散函數(shù)來解決核對最外面電子束縛較弱的問題。一般來說，建議將彌散基組用于處理陰離子的電子親和性、質(zhì)子親和性、反演勢壘以及鍵角等。只對非氫原子加入彌散的s-和p-型高斯函數(shù)，彌散基組用加號“+”表示，如3-21+G；若對重原子和氫原子均加入彌散函數(shù)就用兩個加號“++”表示。1.1.2電子相關(guān)方法由于HF方法為單粒子近似，忽略了電子之間的相互作用，因此會過低的估計鍵長同時會過高地估計振動頻率，因此發(fā)展出了三種考慮電子相關(guān)的方法（或稱post-HF方法）：組態(tài)相互作用（configuration-interaction，CI）M?ller-Plesset(MP)微擾理論密度泛函理論（densityfunctionaltheory，DFT)對典型的分子而言，大部分的基態(tài)能量都源自于單電子HF的能量。MP微擾理論假設(shè)電子相關(guān)的效應(yīng)很小，可以用很小的矯正（即微擾）來得到HF的解。MP方法假設(shè)分子真正的Hamiltonian量可分成兩部分：（其中：代表單電子的貢獻；代表組態(tài)相互作用的貢獻）λ則代表能量和分子波函數(shù)的展開系數(shù)，即：；代入分子的Schr?dinger方程中就可以得到：EHF=E(0)+E(1)+E(2)+…這就是所謂的MP2、MP3、MP4、MP5等等。M?ller-Plesset微擾理論這種方法是將電子相互作用看作是電子密度的函數(shù)，現(xiàn)行的DFT方法通過Kohn-Sham方程將電子能量分成4個部分：E=ET

+EV

+EJ

+EXC，其中ET

代表電子的動能(來自于電子的運動)；EV

代表電子的勢能（包括核與核之間的排斥能、核與電子之間的相互作用能）；EJ

代表電子之間的排斥能；EXC代表電子相關(guān)能。由于對EXC的定義不同便產(chǎn)生了不同的DFT方法，如：LDA、GGA、B3LYP、B3PW91等等。DFT方法可分成三類：局域密度近似（LDA）、廣義梯度近似（GGA）和DFT與HF項結(jié)合的雜化方法：LDA方法只含有與電子密度有關(guān)的項，可以用來預(yù)測塊材料的性質(zhì)，但對于孤立分子的性能預(yù)測效果就不行；GGA方法同時考慮了電子密度和電子密度梯度，用它預(yù)測的鍵離解能與實際僅相差幾個kJ/mol。但GGA方法卻嚴(yán)重地過低估計了一些反應(yīng)的活化位壘，這是因為這種方法忽略了電子的庫侖子相互作用的原因；通過DFT與HF雜化的方法加以解決GGA方法產(chǎn)生的問題。密度泛函理論（DFT）1.1.3閉殼層和開殼層計算大多數(shù)分子具有電子單線態(tài)基態(tài)（電子的總自旋數(shù)為零），這時分子內(nèi)電子完全成對地在限制在軌道內(nèi)，這時就要用閉殼層自旋限制計算。對于那些具有非成對電子的體系，就要采用開殼層非自旋限制方法計算，計算中考慮電子的自旋狀態(tài)（分為自旋向上α電子和自旋向下β電子）。

開殼層非限制性自旋計算適用的體系：具有奇數(shù)個電子的分子(例如一些離子和自由基)；激發(fā)態(tài)；分子有三線態(tài)或更高自旋的基態(tài)(例如O2)；電子對的分離過程(例如分解等)；共振體系，存在電子離域。自旋污染非限制性自旋波函數(shù)并不是自旋算符的本征函數(shù)，因此要對自旋態(tài)進行混合，從而導(dǎo)致自旋污染。在大多數(shù)開殼層系統(tǒng)中，自旋污染是可以忽略的。然而對于未成對電子離域的分子，自旋污染效應(yīng)相當(dāng)明顯，尤其是次高自旋態(tài)。NO的非限制性計算，NO含有一個未成對電子（雙線態(tài)基態(tài)2π），當(dāng)用小基組(STO-3G或3-21G)計算時就呈現(xiàn)出嚴(yán)重的自旋污染。1.2半經(jīng)驗方法半經(jīng)驗方法是對從頭算方法的簡化，如通過限制分子軌道的選擇或僅考慮價電子等再同實驗數(shù)據(jù)進行擬合（如對有機分子的結(jié)構(gòu)和構(gòu)象能量擬合），由此可以大大加快計算的速度。擴展Hückel計算，它忽略了電子之間的相互作用，因而計算速度很快，但精度較差。該方法提供了分子軌道的空間形狀和相對能量的定性估計以及電子密度的空間分布近似。該方法對化學(xué)結(jié)構(gòu)的可視化以及化學(xué)反應(yīng)的前線軌道處理非常有用。忽略微分重疊(Neglectofdifferentialoverlap，NDO)方法，NDO模型忽略了一些但不是全部電子-電子之間的相互作用。Zerner’sINDO方法

(ZINDO/1andZINDO/S)，MichaelZerner’s(Florida大學(xué))開發(fā)的INDO版本，主要用于含過渡金屬的分子體系，ZINDO/1用于預(yù)測分子幾何而ZINDO/S則用于預(yù)測UV光譜。忽略雙原子微分重疊(Neglectofdiatomicdifferentialoverlap,NDDO)方法，該方法基于INDO方法，提供了對完整HF矩陣的最相近的處理。它僅對屬于不同原子的原子軌道使用零微分重疊近似，而保留全部雙中心排斥積分。由于所保留的積分?jǐn)?shù)目太多，因而難以過渡到半經(jīng)驗水平上去。對某些有機小分子和無機絡(luò)合物的非經(jīng)驗NDDO計算，結(jié)果并不令人滿意。改進的忽略雙原子微分重疊(modifiedneglectofdiatomicdifferentialoverlay;MNDO，MNDDO)。它用半經(jīng)驗?zāi)Ｐ吞幚黼p中心排斥積分，其主要特點是把雙原子的電荷分布之間的相互作用看成是兩堆電荷分布的多極距相互作用之和，然后用實驗數(shù)據(jù)優(yōu)化參數(shù)。其參數(shù)僅包含單中心單電子能量和軌道指數(shù)，不僅所用參數(shù)比MINDO/3少，計算結(jié)果也比MINDO/3大有改進，而且機時只增加20％。但它對原子空間排列比較緊密地分子，如四員環(huán)、氫鍵、超鍵化合物、硝基和過氧化物等的預(yù)測則較差。一般地MNDO會過高地估計化學(xué)反應(yīng)的活化位壘。AM1方法(AustinMethod,version1)，對MNDO方法的參數(shù)重新處理，包括核排斥項的改變。盡管AM1方法比MNDO方法精度要高，但它對磷氧鍵、硝基化合物以及過氧化物的鍵的處理不好。PM3、PM5、PM6方法(ParameterisationModel,version3、5和6)，另一種對MNDO方法的參數(shù)重新處理方法，功能類似于AM1，但比AM1方法有明顯的改進。至今還沒有發(fā)現(xiàn)其缺點。DFTB方法（DFT+TB），是將密度泛函方法與緊束縛方法結(jié)合起來，具有處理含有化學(xué)反應(yīng)的較大體系的計算能力。典型的分子力學(xué)方法基本假設(shè)每個原子都可以用一個具有特定質(zhì)量的粒子來代表；化學(xué)鍵可用弦來表示，用以表征兩個原子間相互作用勢能的力常數(shù)，勢能函數(shù)可以描述分子內(nèi)鍵的彎曲、拉伸和扭轉(zhuǎn)以及分子間的一些現(xiàn)象，如經(jīng)典相互作用或范德華(vanderWaals，VDW)力；勢能函數(shù)依賴經(jīng)驗方法如從實驗或從其它計算方法中得到?，F(xiàn)行的分子力學(xué)方法是用一組勢能函數(shù)來表征，用以描述化學(xué)力—力場，這些力場依賴于：原子位移（即鍵長）、原子類型和一些與原子類型及鍵相關(guān)的其它參數(shù)。AMBER(AssistedModelBuildingwithEnergyRefinement)力場,最初用于研究生物分子如蛋白質(zhì)和核苷；CHARMM(ChemistryatHARvardMolecularMechanics)力場，最初也是用于生物分子和藥物研究，但也被用來研究膠束和自組裝大分子；MMx(MM2,MM3,等)

.用以研究非極性小分子的結(jié)構(gòu)優(yōu)化和熱力學(xué)研究MMx力場包含對標(biāo)準(zhǔn)的鍵和鍵角勢能面的二次擬合進行三級和四級校正，因此它可以處理分子振動中的非諧效應(yīng)。不同的MMx版本其參數(shù)也不同，版本越高的MMx其較前一版本而言針對性就越強但適用性就越窄。例如較新的MM4版本就不能夠適用于所有種類的分子；OPLS(OptimisedPotentialsforLiquidSimulations)力場，用于再現(xiàn)生物分子在溶液中的物理性質(zhì)。因為這些力場方法都是針對生物分子的，這里我們不作過多的介紹。但這些方法對含能分子幾何快速的初步優(yōu)化還是有用的。1.4分子動力學(xué)方法(MD)分子動力學(xué)方法是假設(shè)一個體系有N個原子，而體系的狀態(tài)由這N個原子的位置{ri}和動量{pi}或速度{vi}來標(biāo)志，體系的能量為H({ri,pi})，體系的運動方程為：

分子動力學(xué)的主要目的是解上面的方程求得體系狀態(tài)相空間演化的軌跡{ri,pi}t0，{ri,pi}t1，{ri,pi}t2，{ri,pi}t3，……。進而可計算我們感興趣的物理量的值Q({ri,pi})。在實際應(yīng)用中，我們常將上面的哈密頓方程化為下面的牛頓方程，并且用位置ri和速度vi做為描述體系的參量：；式中，V({ri})是原子間相互作用勢，通過解上面的牛頓方程我們可以得到體系在相空間的軌跡，進而求得物理量的平均值[t(1),t(2),t(3),…t(M)]：，（m=1，2，……M）解牛頓方程的方法有很多種，通常有Verlet法則、Verlet速度法則、蛙跳法則、Gear提出的校正預(yù)測法等，最常用的為Verlet速度法則。分子動力學(xué)方法中最關(guān)鍵的是粒子間相互作用。系綜選取根據(jù)研究的體系和方法還可將MD方法分成許多系綜：微正則系綜（NVT、NVE）等壓系綜（NPT、NPH）。分子動力學(xué)計算一般采用周期邊界方法來滿足實際模擬的需要，計算系統(tǒng)中分子間作用力時，采取最近鏡象法，如下圖所示：由于使用了最近鏡象法，在計算時就必須采用截斷半徑(cut-offradius)方法計算遠程的作用力，否則就會因重復(fù)計算粒子所受的力而導(dǎo)致不正確的結(jié)果。在MD模擬中截斷半徑不能大于周期盒邊長的一半。粒子的最近鏡像MD的計算步長選取計算步長的選取，在MD模擬中最重要的工作就是選取合適的步長。選取的標(biāo)準(zhǔn)是既能加快模擬速度又不能失去精度，通常的原則是取小于系統(tǒng)最快運動周期的1/10。例如水分子，分子內(nèi)的運動為鍵長和鍵角的變化，分子間的運動質(zhì)心的運動和分子的轉(zhuǎn)動。分子間的運動源自范德華作用，一般較慢，而分子內(nèi)運動則較快。由紅外譜可知它的最大振動頻率為1.08×1014sec-1，則振動周期為T=1/ν=0.92×10-14sec，那么計算時步長最大可取為約9×10-16sec下表給出了各種模擬體系的計算步長建議：體系運動狀態(tài)步長（sec）簡單原子體系移動10-14剛性分子移動、轉(zhuǎn)動5×10-15軟性分子、限制鍵長移動、轉(zhuǎn)動、扭動2×10-15軟性分子、無鍵長限制移動、轉(zhuǎn)動、扭動、振動10-15或5×10-16MD方法的注意點經(jīng)典分子動力學(xué)方法不能對含有化學(xué)反應(yīng)的體系進行模擬，最近發(fā)展的反應(yīng)力場方法(ReaxFF，A.vanDuin,W.A.Goddard,III)和從頭算的CPMD（CarandParrinelloMD）以及半經(jīng)驗的MD方法可以對含化學(xué)反應(yīng)的體系進行模擬。ReaxFF方法ReaxFF方法是嫁接量子力學(xué)(QM)和經(jīng)典力場(FF)方法的橋梁，它屬于經(jīng)驗方法而更接近于FF方法而不是QM方法。ReaxFF反應(yīng)力場模型可以比QM方法花費低得多的計算機能力來模擬化學(xué)反應(yīng)，其與傳統(tǒng)FF方法相比主要的修改在于：引入了鍵級/鍵距方法，給出了共價鍵級連續(xù)解離的描述；加入了鍵級校正方案，可以允許鍵級方法區(qū)分自由基（離域）和基態(tài)（定域）鍵；在所有的鍵相互作用中都引入了鍵級（包括基于價鍵角和二面角相互作用相）；在所有原子間都進行非鍵相互作用（vanderWaals和Coulomb）計算，因此摒棄了不相容原理；加入了可極化電荷模型。AbinitioMD從頭算分子動力學(xué)（abinitioMD）的計算方法同傳統(tǒng)的MD模擬一樣，都是基于牛頓定律，只不過abinitioMD方法通過解Schr?dinger方程得到原子的力常數(shù)，但abinitioMD方法可以同時得到原子的電子結(jié)構(gòu)，正因為如此，abinitioMD方法非常耗機時?！狟OMD、CPMD半經(jīng)驗分子動力學(xué)的計算方法同從頭算MD模擬一樣，只不過半經(jīng)驗

MD方法中部分解來自于實驗或其它低水平方法，由于半經(jīng)驗MD方法也可以同時得到原子的電子結(jié)構(gòu)，而且需要的資源也遠少于從頭算MD方法，因此在含有化學(xué)反應(yīng)的爆轟模擬上有獨到的優(yōu)勢?！狿M6/MD、DFTB/MD2含能材料的引發(fā)反應(yīng)Coffey:含能材料的引發(fā)反應(yīng)來自于含能材料固體內(nèi)位錯的貫通，當(dāng)沖擊波速度足夠高時，位錯產(chǎn)生的能量會直接傳輸給含能材料分子的內(nèi)在振動模式振動能up-pumping模型指出：沖擊波會產(chǎn)生構(gòu)成晶體的分子的最低振動模式所吸收的激發(fā)聲子熱浴，增加的聲子吸收以及分子間振動能的重新分配會進一步激發(fā)更高的振動模式，最終達到平衡過渡態(tài)，從而產(chǎn)生化學(xué)鍵的斷裂并維持化學(xué)反應(yīng)。值得注意的是，由于爆轟波前緣的傳播速度對化學(xué)反應(yīng)而言太快以至于不能驅(qū)動它們熱激活。Gilman指出，在前緣的壓縮會導(dǎo)致未反應(yīng)物質(zhì)的局部金屬化，這種金屬化由共價鍵的彎曲或由密度達到一定時，它們達到臨界點后形成。在共價性材料中，當(dāng)鍵被拉伸和彎曲時，就會有一個核組態(tài)空間區(qū)域，這時的最高占據(jù)軌道（HOMO）和最低未占軌道（LUMO）的能量差值就會減小?！腆w含能材料會導(dǎo)致金屬化因為爆轟波的前沿會引起未反應(yīng)物質(zhì)的壓縮和剪切，特別是對共價性炸藥材料而言，鍵的彎曲相對而言比鍵的壓縮顯得更為重要，這是因為與之相關(guān)的力常數(shù)會變小。

剪切導(dǎo)致金屬化這一觀點的主要證據(jù)是：通過對一系列不同的氨基-硝基-苯炸藥的HUMO-LUMO差值的研究，其“shake-up”譜與其沖擊感度有直接的關(guān)聯(lián)在體積不變的情況下，鍵的彎曲會引起HUMO-LUMO值變小從而引起在鍵合態(tài)的電子的離域化，這樣會引起快速的化學(xué)反應(yīng)。對TATB而言，在體積接近不變時，臨界剪切（彎曲）角與能隙的消失相對應(yīng)。對非完美含能分子晶體進行的模擬研究，如RDX、PETN等含有孔洞和位錯的情況，都表明當(dāng)RDX晶體中存在單及雙孔洞時，對其進行壓縮它們的光學(xué)帶隙會降低，導(dǎo)致它們的絕緣體—金屬相轉(zhuǎn)變所需的激發(fā)能下降。從直覺上講，非均勻分子固體在猛烈撞擊下會引起原子和電子性質(zhì)的劇烈變化。一個可能的機理就是電子的激發(fā)：電子是屬于輕、快和量子化的粒子，因此它首先對外部的微擾做出響應(yīng)。接下來就會電子激發(fā)，幾步就會引起化學(xué)反應(yīng)的發(fā)生。這包括激發(fā)態(tài)表面的離解、輻射過程以及非輻射能的鈍化。零溫下，能帶閉合也不能爆炸??！

TATB中硝基基團的彎曲振動受剪切的TATB在Fermi面附近的能帶結(jié)構(gòu)β-HMX的點電荷靜電勢3靜高壓下含能材料的行為現(xiàn)在認(rèn)為TATB晶體和硝基甲烷中硝基的彎曲會導(dǎo)致能帶的閉合，考察分子內(nèi)硝基的彎曲在現(xiàn)在實驗所能達到的壓力范圍內(nèi)能帶是否能夠閉合可以通過靜高壓、均勻壓力和非軸向加壓下對硝基甲烷的電子結(jié)構(gòu)計算，得到其光學(xué)帶隙來確定。電子結(jié)構(gòu)的計算可以采用自洽電荷密度泛函—緊束縛理論（SCC—DFTB）得到。這是基于密度泛函理論下對標(biāo)準(zhǔn)緊束縛模型的拓展。自洽部分描述的是總能、原子力和電荷的輸運。實際應(yīng)用中對于每個固定的晶胞參數(shù)集，所有的原子位置均由能量最小化確定，即原子是全弛豫的。SCC—DFTB的緊束縛特點是可以處理大的超晶胞，可包含更多個分子。早期的第一性原理計算則只能計算較小的超晶胞和較低的精度，但都給出了類似的帶隙結(jié)果。以硝基甲烷為例，固定比值連續(xù)降低初始晶胞參數(shù)，壓力P利用低溫公式來估計。對靜水壓而言，當(dāng)體積為原來的50%時，壓力為約50GPa，而HOMO-LUMO減小了約0.6eV，即為原來值的13%。HOMO和LUMO同時增加而其差值卻幾乎隨體積的變化單調(diào)下降（參見下左圖）。在非常高的靜水壓下帶隙的變化值示下面右圖中。壓縮率最高達70%。HOMO-LUMO值的減小，從理論上可以認(rèn)為感度在增加。完美硝基甲烷晶體的原胞的HOMO-LUMO在靜水壓下硝基甲烷HOMO-LUMO值在靜水壓和非靜水壓下隨比容的變化的變化值，突然向下的拐點對應(yīng)高的C-H鍵彎曲C-H鍵形變的機理如下：隨著初始C-N共價單鍵的收縮，電荷被傳輸?shù)紺原子和N原子之間的區(qū)域，隨后C-N鍵變得更強，使得該鍵變?yōu)镃=N共價雙鍵。比鄰的C-H鍵作為電子的供源，至少這些H原子中有一個的部分電子被拉走，導(dǎo)致該鍵變?nèi)酢Ａ硪环矫?，局域在硝基附近的軌道并沒有被壓縮的明顯變形。由Mulliken電荷布居可以看出，H原子總的Mulliken電荷丟失隨每個分子的形變而單調(diào)增加。這些發(fā)現(xiàn)意味著化學(xué)反應(yīng)和共價鍵相離子鍵的轉(zhuǎn)化會導(dǎo)致硝基甲烷的爆轟會始于帶隙的閉合之前。最近的計算研究預(yù)測TATB金屬化的壓力在120GPa，遠高于爆轟壓力。對于硝基甲烷的情況，引起爆轟的原因應(yīng)是質(zhì)子的離解。這就有點意思了，在研究含能材料時，應(yīng)該更注重它在單軸壓縮而不是靜水壓縮下的化學(xué)行為，這與Dick對硝基甲烷的研究結(jié)果相一致。V/V0C1H1aH1bH1cC2H2aH2bH2c0.790.012-0.0100.0160.0160.012-0.0100.0160.0200.600.126-0.016-0.026-0.0260.108-0.019-0.008-0.1550.510.2940.014-0.074-0.0740.168-0.051-0.030-0.245硝基甲烷分子1和2中的C原子和H原子的Mulliken電荷（在單晶中沿b方向壓縮） 4硝胺炸藥的降解動力學(xué)模擬對凝聚炸藥在高密度和高溫下的化學(xué)反應(yīng)機理的詳細(xì)描述對于理解這些材料在燃燒和爆轟條件下反應(yīng)面上所發(fā)生的事件有重要意義。如在沖擊條件下，含能材料會被快速升溫到幾千度并被壓縮到幾百個千巴，體積的減少可達30%。這時會引發(fā)復(fù)雜的化學(xué)反應(yīng)，釋放出的大量能量會導(dǎo)致反應(yīng)持續(xù)下去。很明顯弄清極端條件下各種化學(xué)反應(yīng)機理對于構(gòu)建宏觀水平上連續(xù)爆轟模型有重要意義。靜態(tài)模擬結(jié)果基于自洽密度泛函緊束縛算法（SCC-DFTB）產(chǎn)物H2O、N2、CO2和CO的濃度C(t)隨時間t的演化H2O、N2、CO2和CO的有效反應(yīng)生成速率分別為：0.48、0.08、0.05和0.11ps-1（a）（b）（c）（d）β-HMX撞擊加載下的初始分解動力學(xué)圖1一個分子分解成C3H6N6O6

和CH2N2O2

兩部分，另一個分子分解成C3H6N8O6；CH2和氧原子；圖2顯示C3H6N6O6

進一步分解為C2N4O4H4

，然后又和鄰近的HMX生成的片斷CN2O2H2

反應(yīng)生成C3H6N6O6；圖3顯示生成物中有C2N4O4H4、CH2、

N2O、NO和CH2N2O2

等產(chǎn)物，圖4顯示生成物中有C3H6N6O6、C2N3O2H4、C2N2O2H4、NO2和N2O2；圖5顯示生成物中有C3H5O4H6、N2O2、CH2N2O2、CH2、N2O2和NO；圖6顯示到了反應(yīng)后期生成物基本是一些小分子產(chǎn)物，而且里邊還有一些游離的O原子、H原子和N原子等，這些原子很容易在當(dāng)系統(tǒng)溫度冷卻下來的時候，生成H2,O2,N2,及H2O等小分子量產(chǎn)物。我們在此模擬中尚未發(fā)現(xiàn)有游離碳存在，這可能與我們模擬的時間較短有關(guān)，這個問題將在后續(xù)的工作中進一步研究。問題:

沖擊壓縮模擬在以下方面非常必要:詮釋高溫高壓實驗、炸藥的沖擊感度沿著沖擊路徑研究動力學(xué)和亞穩(wěn)態(tài)解決辦法:新的沖擊模擬技術(shù)(MSST)使得在計算速度上大大加快原子尺度下的量子力學(xué)模擬為這種研究提供了工具得到化學(xué)反應(yīng)的直觀圖像沖擊條件下的從頭算模擬計算原理基于Navier-Stokes方程的NEMD方法結(jié)合沖擊動力學(xué)的三個守恒條件，在沖擊條件下粒子的速度(u)與沖擊波速度(vs)之間的關(guān)系可以寫作：

u=vs(1-ρ0/ρ)，（1）Rayleigh線為沖擊波方向上連接初態(tài)與終態(tài)（Hugoniot）的熱力學(xué)路徑：

p-p0=vs2ρ0(1-ρ0/ρ)，（2）在沖擊條件下遵循能量守恒：e-e0=p0(1/ρ0-1/ρ)+vs2(1-ρ0/ρ)/2，（3）對分子動力學(xué)模擬而言，我們采用拉格朗日量：（4）這里T和V分別代表單位質(zhì)量下的動能和勢能,Q為模擬晶胞尺寸下的質(zhì)量參數(shù)。我們可以發(fā)現(xiàn)當(dāng)粒子速度的一階微分為零時從哈密頓量的形式上方程(3)和(4)是等價的，這是因為T+V=e。由（5）可以看出，當(dāng)粒子速度的二階微分為零時，方程(5)就約化為方程(2)。因此通過方程(5)，我們就可以將沖擊Hugoniot關(guān)系結(jié)合到非平衡分子動力學(xué)中。沖擊模擬可以詮釋與化學(xué)反應(yīng)相關(guān)的熱力學(xué)途徑只能模擬靜態(tài)的熱力學(xué)單態(tài)整個依賴于熱力學(xué)路徑的化學(xué)過程可以得到反應(yīng)只與沖擊波速度和初始條件有關(guān)通過一次沖擊模擬可以獲得大量的熱力學(xué)態(tài)靜態(tài)計算vs.沖擊模擬僅為給定溫度和壓力下的平衡態(tài)根據(jù)產(chǎn)物得到數(shù)據(jù)沖擊波速度在9km/s時的初始分解與產(chǎn)物演化反應(yīng)始于C-H鍵及N-O鍵的首先斷裂；隨著反應(yīng)時間的進行，這些碎片會進一步分解成一個N2O2、C2N2O2H4、NO2

和C2N3O2H4；進而CH2和NO2

從各個碎片中分離出來，分別形成CN2O2H2、C2NH4、CNH2

、N2O、N2O2和CN4O4H2

顯示生成物中有N2、O2、NO2、N2O和游離的O、H、CH、N、C等這些離子在最終會形成H2O、CO、CO2和凝聚碳等。5CO在高溫高壓下的歧化反應(yīng)CO氣體在自然界大量存在，也是炸藥爆轟的主要產(chǎn)物之一。CO作為炸藥的降解產(chǎn)物之一，并且一經(jīng)產(chǎn)生便存在于高溫高壓的環(huán)境之中，在這樣的極端條件下很容易形成更復(fù)雜的化合物，因此能否準(zhǔn)確描述在特定條件下產(chǎn)物的化學(xué)組分以及化學(xué)反應(yīng)過程對于理解炸藥爆轟過程以及爆轟產(chǎn)物狀態(tài)方程具有很重要的意義。CO在沖擊狀態(tài)下的聚合分解反應(yīng)CO在沖擊狀態(tài)下的聚合分解反應(yīng)t=0.16pst=0.18ps

氣體CO在5000K,25GPa條件下的一些產(chǎn)物6極端條件下含能固體HMX結(jié)構(gòu)

與性能的理論研究研究背景：感度分子結(jié)構(gòu)晶體結(jié)構(gòu)晶格聚集態(tài)缺陷空穴起爆機理起爆過程效率設(shè)計HMX研究背景：Bowden[1]和Kanllet[2]等人很早就對感度與分子結(jié)構(gòu)之間的關(guān)系做了大量的工作，認(rèn)為認(rèn)識含能材料的感度不僅要從分子內(nèi)基團的特征和數(shù)目，還必須要深入到分子的電子層結(jié)構(gòu)層次，才能觸及問題的本質(zhì)。一般認(rèn)為C-NO2鍵的斷裂常常被認(rèn)為是硝基化合物分子分解過程中非常關(guān)鍵的一步，但有的硝基化合物分子的分解是首先發(fā)生消環(huán)反應(yīng)（如黑索金(RDX)分子）。因而，確定與C-NO2鍵的強度相關(guān)的因素就具有十分重要的意義。-NO2："致敏"-NH3：

"致鈍"[1]BowdenF.P.,YoffeA.D.,InitiationandGrowthofExplosionsinLiquidsandSolids,CambridgeUniversityPress,Cambridge,1952.[2]KamletM.J.,Proceedings6thSymposium(International)onDetonation,SanDiego,California,ADAO59120,1976.

但是僅從分子角度出發(fā)得出結(jié)論往往并不完全正確，因為不同的晶體構(gòu)型其感度相差很大，比如HMX炸藥有四種不同的晶體結(jié)構(gòu)，其中以β-HMX的感度最低。因此我們要尋求感度和其他因素間的關(guān)系。研究背景：研究背景：

HMX即奧克托金或環(huán)四次甲基四硝基胺(octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine)首先是由W.E.Bachmann和J.Z.Sheehan于20世紀(jì)30年代合成，是生產(chǎn)RDX的副產(chǎn)物，但比RDX有著更高的密度和熔點，自其合成以來以其最優(yōu)的綜合性能、極為廣泛的使用率成為最炙手可熱的猛炸藥，許多含能材料的綜合性能均將其作為參照比較的標(biāo)準(zhǔn)。

-HMX

的晶體及椅子型分子構(gòu)型研究背景：包括2個C4H8N8O8分子，以b為對稱軸，屬于單斜晶系的P21/c準(zhǔn)椅狀的C4H8O8N8－HMX的晶體結(jié)構(gòu)包括8個C4H8N8O8分子，屬于斜方晶系的Fdd2空間群準(zhǔn)船狀的C4H8O8N8研究背景：－HMX的晶體結(jié)構(gòu)包括6個C4H8N8O8分子，以c

軸為對稱軸，屬于P6122空間結(jié)構(gòu)研究背景：－HMX的晶體結(jié)構(gòu)包括兩個C4H8N8O8·0.5H2O分子，具有Pn空間對稱群H2O研究背景：研究背景：在實驗方面：

由于HMX的特殊性，實驗條件上受到很大限制。多集中在對-HMX的高壓壓砧實驗，溫度和壓強對-HMX的光譜振動的影響；還有大量的實驗是對-HMX研究樣品的宏觀性質(zhì)的研究。在理論方面：

理論方面的研究多集中對在-HMX在一定條件下的分子間及分子內(nèi)力場、振動、動能模型的研究。而對HMX四種構(gòu)相的系統(tǒng)研究，尤其是其在極端條件下的系統(tǒng)研究就更少。研究內(nèi)容及方法：

在微觀尺度上對四種結(jié)構(gòu)的HMX在極端條件（高溫高壓）下結(jié)構(gòu)的、電子的、光學(xué)的及熱力學(xué)性質(zhì)進行系統(tǒng)研究，對三種純HMX晶體在高溫高壓下的相變、力學(xué)性質(zhì)以及常溫常壓下-HMX基PBXs在不同晶面上的力學(xué)性質(zhì)進行探討，嘗試在介觀尺度上對添加粘結(jié)劑條件下-HMX的相關(guān)性質(zhì)進行討論。

研究內(nèi)容：研究內(nèi)容及方法：研究方法：密度泛函理論經(jīng)典分子力學(xué)方法耗散粒子動力學(xué)方法使用的軟件：MaterialStudio采用了基于密度泛函理論的平面波贗勢法，使用MS中的Castep模塊；電子之間的交換相互作用嘗試使用了廣義梯度近似（GGA）中的PBE和局域密度近似（LDA）中的CA交換關(guān)聯(lián)勢；離子實和價電子之間的相互作用采用的是Vanderbilt超軟贗勢；布里淵區(qū)積分采用Monkhorst-pack方法，幾何優(yōu)化中采用Broyden-Fletcher-Goldfarb-Shanno(BFGS)

的算法，收斂標(biāo)準(zhǔn)為總能量偏差小于1×10-6eV/atom，原子的最大Hellmann-Feynman力小于0.05eV/?，最大應(yīng)力偏差小于0.01GPa，最大位移偏差小于0.002?。結(jié)果與討論：對-HMX采用4×2×3的K點對布里淵區(qū)進行采樣（網(wǎng)格原點位于Γ點），平面波基函數(shù)的截斷能取500eV；對-HMX采用2×1×4的K點，500eV的截斷能；對-HMX采用3×3×1的K點，500eV的截斷能；對-HMX采用1×2×1的K點，380eV的截斷能。在進行深入的模擬之前，對同一體系采用更高的截斷能和更密的K點采樣進行試驗，結(jié)果表明，這些改變使總能量的變化總小于0.002eV/atom。因此，我們所選取的參數(shù)對于對HMX來說是足夠的。結(jié)果與討論：結(jié)果與討論：β-HMX的體積與壓強的關(guān)系整個范圍內(nèi)，用GGA方法得到的數(shù)據(jù)比實驗值偏高，而用LDA方法得到的數(shù)據(jù)比實驗值偏低，相當(dāng)而言LDA方法的結(jié)果要好一點。低壓下的差別都很大。由于壓強的增加，分子間距不斷減小，分子間相互作用力不斷增強，而使得GGA方法能更好的適合此時的體系，在15到25GPa的范圍內(nèi)GGA的方法比LDA方法的結(jié)果好。27GPa時發(fā)生相變，我們的結(jié)果與實驗結(jié)果不符。結(jié)果與討論：各個方向晶格參數(shù)的變化率和壓強間的關(guān)系，b的壓縮性最大，c的壓縮性最小。表明了在各個方向的可壓縮性有很大的差異。β-HMX各項異性結(jié)果與討論：結(jié)果與討論：N-N鍵的壓縮率最大N-O鍵的壓縮率最小結(jié)果與討論：N7-O3-O4幾乎不隨壓強的變化而變化H13-C2-N4的變化最大結(jié)果與討論：H14-C2-N5-C1的變化最大，17.5%結(jié)果與討論：β-HMX整個布里淵區(qū)價帶和傳導(dǎo)帶的分布都很平坦，這是由-HMX的分子特點決定的。價帶和導(dǎo)帶之間的電子帶隙為3.061eV標(biāo)志著β-HMX是一絕緣體，在40GPa時帶隙值為2.436eV。在0GPa時為傳導(dǎo)帶的G點和價帶的B點的間接帶隙，而壓強在20GPa時為G-G的直接帶隙。隨著壓強的增加，壓強對帶寬的增加的影響越來越小。

結(jié)果與討論：隨壓強的升高帶隙在不斷的減小，使得β-HMX變得越來越不穩(wěn)定，傳導(dǎo)性在不斷加強，表明通過加壓可能使其成為半導(dǎo)體或?qū)w。結(jié)果與討論：0GPa時的部分態(tài)密度。可以判斷β-HMX的基本性質(zhì)：在高價帶，C原子的態(tài)密度值比H、O和N原子的態(tài)密度值要低很多，表明了H、O和N原子在高價帶比C原子的貢獻多很多。在-25到-13eV的能量范圍內(nèi),H、O和N原子存在著明顯的峰值。在傳導(dǎo)帶區(qū)域，峰值主要源于C、H和N原子的p態(tài)。

結(jié)果與討論：和0GPa時相比，25GPa時的態(tài)密度范圍明顯加寬，而且存在很大范圍的負(fù)頻，表明β-HMX晶體的不穩(wěn)定性隨壓強的升高而增強。說明β-HMX隨壓強的增加越來越敏感。

結(jié)果與討論：（a）通過德拜模型方法（b）通過聲子譜方法。所有的焓、熵和溫度乘積、隨溫度的升高而不斷增大，自由能隨溫度升高而不斷減小。結(jié)果與討論：通過兩種方法得到的自由能、焓、熵和溫度乘積、熱容隨溫度的變化趨勢類似，但數(shù)值上相差很大。在文獻[1]中的實驗值是在從310到440K溫度范圍內(nèi)β-HMX的熱容變化為499.43J/mol·K增大到588.37J/mol·K，在相應(yīng)的溫度范圍內(nèi)，用聲子譜法得到的熱容值從589.97J/mol·K增加到775.93J/mol·K，而用德拜模型方法得出的熱容值從1363.59J/mol·K增大到1380.35J/mol·K?？梢耘袛嗤ㄟ^聲子譜方法得到的相應(yīng)的數(shù)據(jù)更可靠。德拜模型的準(zhǔn)諧近似對低對稱性的有機分子晶體可能不太適合。結(jié)果與討論：β-HMX的德拜溫度(a)在0GPa時通過聲子譜方法(b)在0GPa時通過德拜模型方法(c)在0K時通過德拜模型方法。變化趨勢的差距很大，也進一步表現(xiàn)出兩種方法的差異。結(jié)果與討論：在不同壓強下的吸收系數(shù)曲線的形狀類似但對應(yīng)的能量值不同，且隨著壓強的增大發(fā)生了平移。吸收系數(shù)的數(shù)值量級很大(1014nm-1)，表明β-HMX在光的作用下會很不穩(wěn)定且隨著壓強的增加不穩(wěn)定性不斷增強。結(jié)果與討論：結(jié)果與討論：-HMX在0和30GPa時的態(tài)密度結(jié)果與討論：-HMX在0和30GPa時的態(tài)密度結(jié)果與討論：δ-HMX在0和30GPa時的態(tài)密度 對-、-、-HMX的態(tài)密度在費米面附近都有峰值，表明了其對應(yīng)的高價帶的帶結(jié)構(gòu)比較平緩。三種結(jié)構(gòu)HMX中的C原子的態(tài)密度有很大不同，是由它們各自的分子對稱性不同決定的；而苯環(huán)上的N和NO2上的N的態(tài)密度的變化卻不大。在傳導(dǎo)帶范圍內(nèi)，峰值主要由N原子的p態(tài)決定；在能量較高的位置范圍內(nèi)，主要由C和N原子態(tài)密度決定；在費米面附近的總的態(tài)密度的峰值隨壓強的增大而增多。結(jié)果與討論：而-HMX的吸收系數(shù)的峰值隨壓強的升高而減小，且-HMX的吸收系數(shù)在數(shù)值上比-、-HMX的小的多，-、-HMX的吸收系數(shù)峰值為105數(shù)量級，而-HMX的僅為104數(shù)量級，說明在高壓下，-HMX對光的敏感度降低，-HMX比β-、-、-HMX更穩(wěn)定些。-、-HMX的吸收系數(shù)的峰值都隨壓強的升高而增大，對應(yīng)的能量值在減小。結(jié)果與討論：小結(jié)：在低壓下N-N鍵的變化對-HMX穩(wěn)定性的影響最大，而高壓下N-O鍵及O7-N3-N4鍵角對-HMX穩(wěn)定性的影響很大。通過準(zhǔn)諧德拜模型和聲子譜兩種方法研究了-HMX晶體的熱力學(xué)性質(zhì)表明用聲子譜方法得出的結(jié)果比較合理，而準(zhǔn)諧德拜模型方法對低對稱性的-HMX的計算不夠準(zhǔn)確，原因可能是德拜模型采用的近似方法對-HMX這樣的低對稱性的晶體不合適。在0-40GPa內(nèi)并未觀察到HMX的相變，可能是我們的模擬體系不足以在空間群上反映其相變。-HMX對光的敏感程度最小，且隨著壓強的增大光敏感的在不斷下降，而-、-HMX的光感度隨壓強的增大而增強，可以推斷在高壓下-HMX比其他幾種結(jié)構(gòu)相晶體都穩(wěn)定。

在MS程序中的Discover模塊下的恒溫恒壓（NPT）的分子動力學(xué)方法模擬-HMX、-HMX和-HMX在不同條件下進行研究，采用COMPASS(condensed-phaseoptimizedmolecularpotentialsforatomisticsimulationstudies)

力場，采用NPT系綜中的Andersen的調(diào)溫方法及Berendsen調(diào)壓方法。

三種結(jié)構(gòu)的HMX的初始結(jié)構(gòu)如前。模擬的超晶胞包括96個HMX分子，2688個原子，對-HMX、-HMX和-HMX相應(yīng)的盒子大小分別為由4×4×3、2×2×3、4×2×2的超晶胞組成，保證分子間作用勢的截斷半徑在10.5?。結(jié)果與討論：HMX的溫度和壓力相變從左到右依次為-HMX、-HMX和-HMX晶體的超晶胞結(jié)構(gòu)示意圖

模擬條件為：-HMX在298K時0到31GPa及0GPa時溫度從5到400K；-HMX在298K時壓強從22到40GPa及0GPa時溫度從350到480K；-HMX在298K時壓強從20到40GPa及0GPa時溫度從430到555K。所有模擬均采用2×105時間步，每個時間步為1fs，其中前1×105時間步用來平衡結(jié)構(gòu)，后1×105時間步用來收集結(jié)果。在平衡結(jié)構(gòu)時，為了保證HMX晶體的溫度和壓強能與模擬外部的溫度及壓強保持一致，每間隔1000時間步標(biāo)定一次。相互作用是由原子核在截斷范圍內(nèi)的最鄰近原子決定的。所有的計算考慮周期性邊界條件，非鍵相互作用分別采用Ewald和基于原子的方法來控制的庫侖力和長程范德瓦爾斯力。結(jié)果與討論：a/?b/?c/?α(o)β(o)(o)V(?3)(g/cm3)β-HMXExpt.[20]6.5411.058.7090124.390259.69351.89369Thiswork6.5811.128.7690124.390264.91351.85638α-HMXExpt.[5]15.1423.895.91909090267.20251.84048Thiswork15.3124.175.98909090276.60751.77790δ-HMXExpt.[8]7.717.7132.559090