NTO結晶形貌的預測

馬 松，袁俊明，劉玉存，常雙君，王建華，于雁武

（中北大學化工與環(huán)境學院，山西 太原 030051）

引 言

3－硝基－1，2，4－三唑－5酮（NTO）是一種性能優(yōu)良的高能鈍感炸藥［1］，其爆轟性能和感度大小與其晶體形貌有一定的關系。NTO分子具有偶極子的特性，其結構的共價鍵具有方向性和飽和性，生長基元在結晶方位和晶體各個面族上附著的穩(wěn)定性影響到結晶形貌［2］，且NTO晶貌特征還與生長時的物理、化學條件有關，因此研究NTO晶體形貌需要把晶體內(nèi)部結構同生長條件密切結合起來。目前關于晶體生長理論和晶貌預測的研究很多，如Singh M R［3］等人根據(jù)鄰苯二甲酸氫鉀其晶體可能的生長面及其面間距，用群體平衡模型（M－PBM）推出 KHP的晶體形貌；WANG Yan－chao［4］等人基于粒子群優(yōu)化算法，計算出單晶硅具有低自由能的生長面，并根據(jù)其對稱性信息和各晶面自由能預測了硅的生長晶習；Ulrich J［5］等人通過修正附著能法預測了苯甲酸在水中的重要生長面、重結晶形貌，其結果得到了實驗驗證；李鑫［6］等人結合實驗與模擬計算，研究了添加劑對重結晶HMX晶體形貌控制的影響。

為了研究NTO晶貌形成的內(nèi)在機理、微觀尺度上溶劑分子與NTO晶面的作用機制和對其結晶形貌的影響，首先用附著能（AE）模型、分子動力學（MD）方法對NTO在真空、水及甲醇中的吸附能和晶體形貌進行了研究，并從NTO分子的極性、與溶劑所成氫鍵的角度分析了控制晶體形貌的內(nèi)在和外在因素，為NTO重結晶工藝中溶劑的選擇提供參考。

1 真空中NTO晶貌的預測

1.1 AE晶體生長模型

晶體生長的附著能（AE）模型是在周期性鍵鏈（PBC）模型的基礎上發(fā)展起來的，PBC模型假設晶體表面成鍵所需時間與鍵合能成反比，晶面生長速度與鍵合能成正比，晶體結構由PBC組成，且晶體生長最快的方向是化學鍵最強的方向。而AE模型在此基礎上定義了晶片能Ea為生長出一層厚度為dhkl的晶片釋放出的能量，附著能Eatt為晶片附著在晶體表面（hkl）所釋放出的能量，兩者之和為晶體的晶格能Elatt［6－9］。

對于有機分子，晶格能可通過非成鍵原子間作用力函數(shù)和分子間作用綜合計算求得，對于晶片能和附著能，可由中心分子與晶片內(nèi)、外的所有分子總作用力計算得出。

Berkovitch－Yellin［10］對上述方法進行了改進，指出各晶面的相對生長速率（Rij）正比于附著能Eatt。

最低附著能的晶面生長得最慢，所以具有最大的形態(tài)學特征，從而決定了晶體的微觀形貌。該方法是基于晶體對稱性、分子間鍵鏈特性和晶面附著能來計算、推出有機分子在真空和有機溶劑中的附著能、晶貌特征。

1.2 NTO晶體模型的構建

在CCDC劍橋晶體結構數(shù)據(jù)庫中可查出計算所需的NTO晶胞參數(shù)，共有α和β兩種晶型，因后者密度較高，在起爆藥內(nèi)應用較多，故以后者為研究對象，β型NTO晶體為單斜晶系，屬于CS點群，其晶胞參數(shù)為：a＝9.322，b＝5.4791，c＝9.0685，α＝γ＝90°，β＝101.322°。

根據(jù)β型NTO晶體中原子分數(shù)坐標和對稱性信息在material studio 6.0中建立原胞，并在Discover模塊下進行5000步的幾何優(yōu)化，用優(yōu)化后的原胞建立3×3×3的超晶胞，為進一步消除NTO超胞中不合理的能量和構象，用Forceit模塊下退火（Annealing）操作對其進行300～500K升、降溫循環(huán)處理，此過程總步數(shù)為10 000步，每隔1 000步保存一次構象，采用Compass力場，NPT系綜，Andersen控溫方法，各原子起始速度按Maxwell分布取樣，Velocity Verlet算法標定原子速度，分別用Atom based和Ewald方法計算范德華和靜電作用，原子截斷半徑取9.5?，緩沖寬度取0.05［11］。經(jīng)此過程可得到結構合理、勢能較低的構象文件，選取最低勢能構象文件作后續(xù)計算，圖1為優(yōu)化后的NTO超晶胞。

圖1 優(yōu)化后的NTO超晶胞Fig.1 The optimized supercell of NTO

1.3 真空中NTO的晶貌

預測NTO真空中晶貌的主要目的是從晶體的周期性鍵鏈（即分子間作用）和自身附著能的角度分析無溶劑影響下晶貌形成的內(nèi)部因素，在此基礎之上用晶面與溶劑分子的吸附能修正NTO晶體重要生長面的附著能，以此預測和研究NTO在水及甲醇中的晶貌。

在Morphology模塊下，用1.2節(jié)中的力場類型、熱力學系綜、截斷半徑等控制方法，對所得豫弛后的構象用附著能法預測NTO晶體形貌，結果見圖2（a），用Crystal Graph法計算優(yōu)化后NTO原胞內(nèi)分子基團間的相互作用，見圖2（b），虛線為分子間力，強弱為藍＞紫＞紅（藍線方向的晶面生長最快，已消失），NTO重要晶面的生長信息見表1。

圖2 真空中NTO晶體形貌（a）和分子間作用力（b）Fig.2 Morphology（a）and molecular interaction（b）of NTO in vacuum

表1 真空中NTO重要晶面的生長信息Table 1 Growth of important surface of NTO in vacuum

從圖2可見，NTO分子在（010）方向形成了一種偶極生長基元，－NO2基團電子云密度大，為負極面，與三唑環(huán)上的正電荷中心（距離C原子較近）有很強的分子間力，見表1中的Eatt（Elec），（010）晶面的法線與－NO2基團伸展方向夾角小，受此極化作用影響，分子在（010）面附著效率較高，故（010）面的生長速率較（100）和（001）面快，且晶面面積較小，相對于其他快速生長、或已消失的面仍頑強顯露。

（001）面的法線幾乎垂直于－NO2基團和三唑環(huán)所在平面，分子的附著主要靠共軛體系間微弱的范德華力作用，晶面生長速率較（010）面低且顯露面積大；（100）面法線亦與－NO2基團伸展方向夾角大且平行于三唑環(huán)共軛平面，另外NTO分子在此方向上側向排布，附著能主要由羰基與三唑基團間的親和作用提供，致使分子的附著效率最低，其晶面所占面積最大，成為NTO形態(tài)學上最重要的生長面。

相對于已消失、或顯露面積小的生長面，（010）、（100）和（001）晶面為重要生長面，因Rij∝Eatt，故真空中各晶面生長速率（010）?（100）≈（001），其形貌為長四棱柱。據(jù)分析可知，真空中NTO晶體形貌主要受其分子環(huán)上鍵鏈作用弱的（100）和（001）晶面決定，其次是極性較強的（010）晶面。

2 溶劑中NTO晶貌的預測

2.1 吸附模型的構建與優(yōu)化

溶劑分子對NTO晶貌的影響是通過與其晶面的吸附作用實現(xiàn)的，故需建立水、甲醇與NTO晶面的吸附模型，再求解其吸附能。水和甲醇分子的模型由Sketch建立后，也采用1.2節(jié)中方法進行幾何構型和能量的豫弛。

由表1可知，（010）、（001）和（100）晶面為三組重要生長面，其面積占晶體總表面積的92.78%，為減小計算量，只計算這三組晶面對水和甲醇分子的吸附能。

對退火所得的3×3×3超晶胞（010）、（001）和（100）面進行切割，切割厚度為一個分子層，再對此二維結構加20?的真空層，由此形成了具有吸附面的層晶。分別往各個吸附面加入水和甲醇分子，用文獻［12］中的方法調(diào)整合理的構型，以形成6個初始吸附模型，對每種模型進行幾何優(yōu)化，方法與前述相同。經(jīng)優(yōu)化后的水、甲醇與NTO各晶面的吸附模型見圖3（虛線為氫鍵），其中（100－H2O）吸附模型與文獻［2］中的NTO與水全優(yōu)化構型相符。

圖3 NTO晶面與水分子（a）和甲醇分子（b）的吸附模型Fig.3 Sorption models of NTO crystal surfaces with water（a）and methanol（b）

2.2 吸附過程的動力學計算

固定除表層以外的分子層以減少計算量（溶劑分子與固定層距離大于9.5?），用Forceit模塊下的動力學過程對優(yōu)化的6種吸附模型進行NPT系綜下的MD計算，總時間為150ps，步長為1fs，其他控制參量與退火過程相同，MD過程結束后，用各模型能量穩(wěn)定后的構型（后50幀）分別對整個體系、溶劑分子和晶面的能量求均值，再用吸附能計算公式（3）求解吸附能［6－9］，計算結果見表2。

式中：Etotal為吸附體系的總能量（kJ／mol）；Esolv為溶劑分子單點能（kJ／mol）；Esurf為晶面的表面能（kJ／mol）。

表2 各晶面與溶劑分子的吸附能Table 2 The sorption energy between surface and solvent

2.3 溶劑中NTO的晶貌

為研究溶劑及其他因素對NTO晶體形貌的影響，需要用表2中的吸附能對表1中的附著能進行修正（本研究僅考慮溶劑因素）。因各溶劑對晶體不同NTO晶面的附著能影響不同，可用重要生長面（010）、（001）和（100）對溶劑分子（水和甲醇）的吸附能對此三組晶面上NTO分子的附著能予以修正，計算公式為［8－9］：

式中：Eatt′為修正后的附著能（kJ／mol）；Ebind為晶面對溶劑分子的吸附能（kJ／mol）；Eatt為真空中的附著能（kJ／mol）。另外由AE模型理論可推斷出修正后的生長速率、面心距和附著能亦符合R′ij∝D′ij∝E′att，三者的計算結果見表3，可見溶劑分子對晶貌的影響是通過與晶面的吸附來影響NTO分子在各面的附著能量與效率，溶劑與面相吸附能愈大，面生長愈慢，反之亦然。

表3 水和甲醇中NTO主要晶面的附著能Table 3 Attachment energy of important surface of NTO in water and methanol

根據(jù)水、甲醇與不同晶面吸附能，修正各晶面的附著能值，為各晶面指定相對生長速率和中心到面的距離，再通過Wulff圖對以上信息進行分析，推導出晶體形態(tài)，并與文獻［13］中NTO在水和甲醇中的結晶形態(tài)電鏡掃描圖進行對比，結果見圖4和圖5。

圖4 NTO 在水中的晶貌預測（a）和SEM 圖（b）［13］Fig.4 Prediction of morphology on NTO in water（a）and SEM image（b）

圖5 NTO在甲醇中的晶貌預測（a）和SEM 圖（b）［13］Fig.5 Prediction of morphology on NTO in methanol（a）and SEM image（b）

據(jù)表1～表3可知，NTO各晶面會通過對水和甲醇分子的吸附而降低自身的附著能。文獻［6，8－9］亦顯示，極性溶劑與有機分子的相互作用中，氫鍵的成鍵強度較范德華力和靜電力顯著。

綜合圖4和表3可知，（010）晶面上的硝基分布密集，水分子中的氫易于硝基氧締合成氫鍵，且成鍵密度較高，所以該面對水分子吸附能最大，即水分子對（010）面附著能的抑制最明顯；水分子中的氫、氧原子能與NTO分子上的9位氧、11位氫原子［2］間形成兩個氫鍵且方向垂直于（100）面，故與（100）面吸附能也較強，但該面在真空中的附著能基數(shù)較小，水的抑制作用較（010）面次之；（001）面與水的吸附能由氫鍵和范德華力共同作用（該面上氫、氧原子分布少，氫鍵對吸附能貢獻?。?，且量值最小，水分子對附著能抑制作用最小。故水中各面生長速率為（010）＞（001）＞（100），其結晶形態(tài)為扁平的六面體。

從圖5和表3可見，甲醇分子同各面的吸附機理類似于水，但由于甲基推電子集團使羥基周圍電子云密度高，分子偶極距大［7－9］，與各面形成的氫鍵作用較強且鍵能差距較小，致使（010）和（001）面的附著能降低的量值相近；另外在（100）面所成的氫鍵無水分子與該面所成氫鍵的特殊性，故對（100）面的抑制強度與水相當。故甲醇中各面的生長速率（010）＞（100）≈（001），其結晶形態(tài)為四棱柱。

總之，溶劑分子對晶貌的影響由NTO分子內(nèi)氧原子的排布、數(shù)目，溶劑分子偶極距、羥基數(shù)目所決定。用修正后的AE模型所預測的晶貌與文獻［13］（P－84）中掃描電鏡所觀測的實際晶貌能較好地匹配，說明預測是較成功的。

3 結 論

（1）利用分子動力學和修正附著能的方法可對NTO晶體在真空、水和甲醇環(huán)境中的晶貌進行預測，其結果和文獻中的掃描電鏡圖能較好地匹配。

（2）真空中，NTO晶體的主要生長面為（010）、（100）和（001）面，因附著能隨晶面法向極性的增強而升高，故非極性面（100）和（001）生長最慢，有硝基顯露的極性面（010）次之。

（3）水和甲醇通過與NTO分子環(huán)上的硝基或羰基氧形成的氫鍵而降低各晶面附著能，且硝基與醇羥基所成氫鍵最強。表現(xiàn)為水中各面生長速率（010）＞（001）＞（100），甲醇中（010）＞（100）≈（001），且總體上甲醇中的生長速率稍小。

［1］ Nandi A K，Singh S K，Kunjir G M，et al.Assay of insensitive high explosive 3－nitro－1，2，4－triazol－one（NTO）by acid－base titration ［J］.Central European Journal of Energetic Materials，2013，10 （1）：113－122.

［2］ 方國勇，徐麗娜，肖鶴鳴，等.3－硝基－1，2，4－三唑－5－酮與NH3及H2O分子間相互作用的理論研究 ［J］.化學學報，2005，63：1056.FANG Guo－yong，XU Li－na，XIAO He－ming，et al.Theoretical study on intermolecular interactions of 3－nitro－1，2，4－triazol－5－one with NH3and H2O ［J］.Acta Chimica Sinica，2005，63：1056.

［3］ Singh M R，Ramkrishna D.A comprehensive approach to predicting crystal morphology distributions with population balances ［J］.Cryst，Growth Des，2013，13：1397－1411.

［4］ WANG Yan－chao，LüJian，ZHU Li，et al.Crystal structure prediction via particle－swarm optimization［J］.Physical Review B，2010，82：094116.

［5］ Ulrich J.Predicting crystal morphology grown from solution ［J］.Chem， Eng， Technol，2012，35 （6）：1009－1012.

［6］ 李鑫，陳樹森，李麗潔，等.添加劑對HMX重結晶晶體形貌的影響 ［J］.火炸藥學報，2011，34（3）：15－19.LI Xin，CHEN Shu－sen，LI Li－jie，et al.Influences of additives on HMX crystal morphology ［J］.Chinese Journal of Explosives and Propellants，2011，34（3）：15－19.

［7］ ZHOU Kun，LI Jun，LUO Jian－h(huán)ong，et al.Crystal growth，structure and morphology of rifapentine methanol solvent ［J］.Chinese Journal of Chemical Engineering，2012，20（3）：602－607.

［8］ CHEN Jian－xin，WANG Jing－kang，ZHANG Ying，et al.Crystal growth，structure and morphology of hydrocortisone methanol solvate［J］.Journal of Crystal Growth，2004：265－273.

［9］ 陳建新，王靜康，尹秋響，等.氫化可的松形貌預測［J］.天津大學學報，2006，39（增刊）：3－6.CHEN Jian－xin，WANG Jing－kang，YIN Qiu－xiang，et al.crystal morphology prediction of hydrocortisone.［J］.Journal of Tianjin University，2006，39（Suppl）：3－6.

［10］Berkovitch－Yellin Z.Toward an ab initio derivation of crystal morphology［J］.J Am Chem，Soc，1985：107.

［11］趙樹森，付一政，梁曉艷，等.HTPB／增塑劑共混物的介觀動力學模擬 ［J］.高分子材料科學與工程.2011，27（5）：186－187.ZHAO Shu－sen，F(xiàn)U Yi－zheng，LIANG Xiao－yan.et al.Mesoscale dynamic simulation on HTPB／plasticizers blends［J］.Polymer Materials Science and Engineering，2011，27（5）：186－187.

［12］王志明.羥基磷灰石晶習及降解特性的研究 ［D］.天津：天津大學，2007.

［13］Teipel U.Energetic Materials［M］.Weinheim：Wiley－Vch Verlag GmbH ＆Co.KGaA，2005：84－139.