周 濤，陳 芳，李 軍，曹端林，王建龍

（中北大學化學工程與技術學院，山西 太原 030051）

1 引言

5，5′?聯四唑?1，1′?二氧二羥銨（TKX?50）是一種四唑類含能離子鹽，區(qū)別于傳統(tǒng)的含能材料，TKX?50分子結構中并沒有硝基致爆基團，但是爆速和爆壓等爆轟性能均超過了奧克托今（HMX）［1］，主要原因是環(huán)形骨架中存在高能的N—N、N=N 鍵。自從TKX?50被合成后，國內外針對其衍生物設計以及合成工藝路線優(yōu)化進行了廣泛的研究［2-6］。一般來說，含能材料的晶體形貌是影響其性能的重要因素之一［7］，在實際應用中經常期望得到特定的形貌，以滿足其安全性能需求。例如，晶體形貌規(guī)則且接近球形，可以提高炸藥的裝藥密度和流散性，進而改善藥柱的力學性能和機械感度。徐容等［8］利用溶劑?非溶劑重結晶法得到粒度跨 度 小、撞 擊 感 度 低 的2，6?二 氨 基?3，5?二 硝 基 吡嗪?1?氧化物（LLM?105）晶體；周誠等［9］在二甲基亞砜/水（DMSO/H2O）混合溶劑中重結晶得到表面光滑、熱感度較低的1，1?二氨基?2，2?二硝基乙烯（FOX?7）晶體；許誠等［10］采用降溫重結晶法獲得粒度跨度小、表面光滑的TKX?50 晶體，降低了其機械感度。

近年來，計算機科學的蓬勃發(fā)展，使得人們逐漸掌握了在原子分子水平上模擬晶體生長過程的手段。其中，分子動力學（MD）方法已經廣泛地被應用于研究炸藥晶體生長過程中的溶劑效應。Duan 等［11］借助MD 方法解釋了丙酮溶劑對HMX 晶體形貌的影響。Chen 等［12］通過MD 方法計算了黑索今（RDX）晶體在丙酮溶劑中的生長形貌；Lan 等［13］采用MD 方法模擬了外部生長環(huán)境對六硝基六氮雜異伍茲烷（HNIW）晶體形貌的影響，外部生長環(huán)境包括溶劑、溫度和過飽和度等。關于TKX?50 在一元溶劑中的晶形預測已有報道［14］，但是TKX?50 在二元混合溶劑中的晶體形貌尚未被系統(tǒng)研究，混合溶劑體積比與溫度對TKX?50 晶體形貌的影響情況尚不清楚。而TKX?50 在大多數常用溶劑中的溶解度非常差，如二甲基甲酰胺（DMF），甲醇，乙醇，甲苯，乙酸乙酯等，但在甲酸和水中具有更好的溶解性［15］。因此，本研究采用MD 方法模擬TKX?50在甲酸/水混合溶劑中的晶體形貌，采用附著能（AE）模型預測TKX?50 在真空中的晶體形貌，預測TKX?50 在不同溫度、不同體積比的甲酸/水中的結晶形貌，并通過實驗結果來驗證模擬的準確性；最后進行徑向分布函數分析，以期了解TKX?50 晶面與混合溶劑之間存在的相互作用力。

2 計算細節(jié)

2.1 AE 模型及其修正

AE 模型是在周期性鍵鏈（PBC）理論基礎上建立起來的模型［16-17］，它根據晶體對稱性和分子間鍵鏈性質計算晶體的附著能。該理論認為各個晶面相對生長速率與晶面附著能的絕對值成正比［18］：

式中，Rhkl為相對生長速率；Eatt為晶面的附著能，kJ·mol-1。

晶面附著能絕對值越高，晶面生長速率越快，晶面生長趨于減小或消失；晶面附著能絕對值越低，晶面生長速率越慢，在最終的晶體形態(tài)中越容易得到顯露。通過評估晶體不同晶面的相對生長速率Rhkl，可以應用AE 模型預測晶體的習性。然而，該模型的準確度有待商榷，因為AE 模型忽略了外部結晶條件，并不能準確反映溶液結晶過程中晶體的實際生長過程。因此，對AE 模型進行修正，修正后的附著能可以通過下式來計算［11，13，19-20］：

式中，E代表修正后的附著能，kJ·mol-1；S用來描述表面的特征，其定義為：

式中，Aacc為單位晶胞（h k l）面的溶劑可及面積，?2；Ahkl為單位晶胞（h k l）面的面積，?2。

ES代 表 溶 劑 對 晶 面 的 影 響 能 力，kJ·mol?1，可 以 通過下式來計算：

式中，Abox是超晶面模型的面積，?2；Eint是溶劑?表面的相互作用能，kJ·mol-1，其計算公式為：

式 中，Etot是 溶 劑 層 和 晶 面 層 的 總 能 量，kJ·mol-1；Esur（Esol）是去除溶劑層（晶面層）的晶面層（溶劑層）的總能量，kJ·mol-1。

2.2 模擬過程

TKX?50 的初始晶胞結構來源于單晶衍射數據［1］，其分子與晶胞結構如圖1 所示。TKX?50 屬于單斜晶系，P21/C 空間群，晶胞參數為a=5.44 ?，b=11.75 ?，c=6.56 ?，α=γ=90°并且β=95.07°，每個原胞中含有2 個TKX?50 分 子。力 場 選 擇PCFF（Polymer Consis?tent Force Field），其參數得到了驗證［14，21-22］，適用于TKX?50 的模擬研究。在模擬過程中，非鍵相互作用，包括范德華力［23］和靜電力［24］，分別通過Atom?based和Ewald 求和方法計算。Atom?based 方法的截斷半徑設定為12.5 ?，Ewald方法的精確度為0.004 kJ·mol-1。

圖1 TKX?50 的分子和晶胞結構圖（C，H，O，N 分別由灰色，白色，紅色和藍色表示）Fig.1 Molecular and unit cell structures of TKX?50（C，H，O，N are represented by gray，white ，red，and blue colors，respectively）

使用AE 模型，預測真空中TKX?50 的晶體形貌，確定形態(tài)學上重要的晶面。然后對TKX?50 的重要晶面進行切割，將其構造為3×3×3 的周期性超晶胞結構模型，并對超晶胞結構模型進行優(yōu)化。隨后，構建溶劑層模型，溶劑體系為不同體積比的甲酸/水（1/4，1/3，1/2，1/1 和2/1）混合溶劑。對溶劑層進行幾何優(yōu)化，將優(yōu)化后的溶劑層沿著c軸與晶面層對接，構建晶體?溶劑雙層模型，以此來探究溶劑對晶體形貌的影響。晶面層與溶劑層之間的真空距離為3 ?，而溶劑層上方空出50 ? 的真空距離，以便消除上下晶體表面的影響。對雙層模型進行20000 次的迭代優(yōu)化，然后進行MD 模擬，MD 模擬期間要固定晶面層。MD 模擬總時間設定為300.0 ps（300000 fs），時間步長為1 fs，系綜選擇NVT，溫度由Andersen 恒溫器控制［25］。當溫度和能量的波動小于10%時，整個體系可以認為達到了平衡。

3 結果與討論

3.1 真空中TKX?50 的結晶形貌

TKX?50 晶體單晶衍射數據與優(yōu)化后的晶胞參數列于表1。從表1 看到，只有c軸長度偏差較大（6.95%），但仍在可接受的偏差范圍內，這說明本次模擬力場選擇PCFF 是合適的。圖2 為TKX?50 分子之間的相互作用。橢球體為TKX?50 晶體的生長基元，橢球體中心代表生長基元的質心。TKX?50 晶體從生長基元開始向外生長，從質心向外擴散的線是TKX?50 分子之間的相互作用，同時也是晶體生長的驅動能量，其能量大小與長度之間的關系如表2 所示。藍色線代表能量較大的強鍵，其值為-460.62 kJ·mol-1；紅色線代表能量較小的弱鍵，其值為-163.42 kJ·mol-1。晶體在真空中的生長形貌，是由其表面鍵能決定的，表面鍵能越大，生長速度越快，因此TKX?50 分子沿著藍色線方向比紅色線方向生長速率快。生長速率較快的面趨于消失，而生長速率較慢的面容易保留下來，最終形成TKX?50 在真空的晶體形貌。

表1 TKX?50 晶胞參數的實驗數據和優(yōu)化數據的對比Table 1 Comparisons of experimental and optimized cell pa?rameters for TKX?50

圖2 TKX?50 分子之間的相互作用Fig.2 Interactions between TKX?50 molecules

表2 TKX?50 分子間相互作用大小Table 2 Intermolecular interaction sizes of TKX?50

任曉婷等［21］預測TKX?50 真空中晶體晶形為不規(guī)則長塊形，主要包括4 個晶面，分別是（0 2 0）、（1 0 0）、（0 1 1）和（1 1 0）晶面。Xiong 等［22］預測TKX?50 真空中晶體晶形為片狀，主要晶面為（0 2 0）、（1 0 0）、（0 1 1）、（1 1 0）、（1 1 1ˉ）和（1 2 1ˉ）。劉英哲等［14］預測TKX?50真空中晶體晶形呈長片狀，由5 個生長晶面構成，即（0 2 0）、（1 0 0）、（0 1 1）、（1 1 0）和（1 1 1ˉ）。本模擬預測的TKX?50 在真空中的形貌如圖3 所示，主要有5 個重要生長面，分別是（0 2 0）、（1 0 0）、（0 1 1）、（1 1 0）和（1 1 1ˉ），晶體晶形為不規(guī)則長塊形。本研究模擬的結果與文獻［14］中的模擬結果相一致。縱橫比定義為晶體習性的最長和最短直徑之間的比值?？v橫比越接近1，預測的晶體形貌越接近球形。本研究模擬預測TKX?50 真空形貌的縱橫比為1.98，其形貌參數列于表3。其中，（0 1 1）面的面積比最大，達到了35.65%，具有最大的形態(tài)重要性。（1 1 1ˉ）面的附著能為-312.63 kJ·mol-1，其絕對值相對最大，而且（1 1 1ˉ）面占總面積比僅有0.36%。因此，（1 1 1ˉ）面的生長速率快于其它重要面，趨于消失。

圖3 TKX?50 在真空中的形貌預測圖Fig.3 Morphology prediction diagram of TKX?50 in vacuum

表3 真空中TKX?50 形態(tài)重要晶面的參數Table 3 The parameters of morphologically important faces of TKX?50 in vacuum

TKX?50 分子的羥銨陽離子和聯四唑陰離子通過強靜電作用結合，不同的分子排布方式可能會導致不同的電荷分布，最終影響晶面與溶劑之間的相互作用。TKX?50 各個重要晶面的分子堆積結構見圖4?？梢钥吹?，在（0 2 0）晶面上有排列規(guī)則的聯四唑陰離子顯露。（1 0 0）晶面上聯四唑陰離子環(huán)對立分布，空間位阻較大。（1 1 0）和（1 1 1ˉ）晶面的分子堆疊相對分散，既有聯四唑陰離子顯露，也有少量羥銨陽離子顯露，陽離子顯露部位空間位阻較小。因此溶液中的溶質分子能夠很容易地吸附到（1 1 0）和（1 1 1ˉ）晶面上，導致晶面快速生長。晶面的粗糙程度會影響晶面的生長速率。粗糙晶面具有較多的生長臺階和扭結點位，容易吸附溶液中的溶質分子［26］，因此具有較快的生長速率。參數S可以用來描述晶面的粗糙程度，計算方法見公式（3）。表4 列出了TKX?50 不同晶面的S值。（1 1 0）和（1 1 1ˉ）晶面的S值分別為1.45 和1.35，相對其它晶面的S值更大，說明（1 1 0）和（1 1 1ˉ）晶面相對其它晶面更粗糙，在真空中生長更快。（1 0 0）晶面的S值最小，表明（1 0 0）晶面相對平坦，在真空中生長較慢。此外，由于（1 1 0）晶面的S值最大，具有最粗糙的表面特征，可以預測，（1 1 0）表面吸附點位較多，也可能會與溶劑分子產生較強的吸附作用。

圖4 TKX?50 重要晶面的分子堆積結構Fig.4 Molecular stacking structure of TKX?50 habit faces

表4 TKX?50 重要生長面的溶劑可及面積和表面面積Table 4 Solvent?accessible areas and surface areas of mor?phologically important crystal faces for TKX?50

從前面的分析來看，（1 1 0）晶面由于其表面的結構特征，可以對溶質分子和溶劑分子產生較強的吸附能力。但是，晶面更容易吸附溶質分子還是溶劑分子，僅分析表面結構是不夠的。溶液結晶是一個復雜的相轉變過程。一方面，溶劑吸附于晶體降低了界面能，使晶面由光滑面轉變?yōu)榇植诿妫@會促進晶面的生長；另一方面，溶劑的優(yōu)先吸附占據了晶面的生長活性位，溶質生長必須先克服溶劑的脫附能壘，這會阻礙晶面的生長［27］。溶質分子和溶劑分子在界面上競爭吸附，其吸附強弱能力需要從多方面來分析，比如晶面結構，溶質和溶劑分子結構，表面靜電勢等。本課題組前期工作 分 析 了TKX?50 晶 面 的 表 面 靜 電 勢［28］，認 為 除 了（1 0 0）晶面，其它晶面都為極性面。晶面的正電荷比較多，因此溶質中的陰離子相較于陽離子來說更容易吸附到晶面上。此外，選擇具有負電子基團或富電子芳環(huán)的溶劑分子，更容易控制TKX?50 晶面的生長速率，從而控制TKX?50 的晶體形貌。

3.2 不同體積比的甲酸/水混合溶劑對TKX?50 晶體形貌的影響

溶質分子從溶液主體擴散到晶體表面并且長入晶體，促進了晶體的生長。然而，溶劑會優(yōu)先占據晶面的生長活性位點，溶質分子吸附到晶面必須克服溶劑的脫附能壘［29-30］。因此，溶劑在晶面上的吸附會影響晶面的生長速率，促使不同晶面的各向異性生長，導致晶體的生長形貌不同。以（0 1 1）和（1 1 1ˉ）晶面為例，MD 模擬（0 1 1）和（1 1 1ˉ）晶面上甲酸/水的分子分布如圖5 所示?？梢钥吹剑? 1 1）晶面相對平坦，具有較大的空間位阻，這使得甲酸和水分子難以與（0 1 1）晶面上的羥銨陽離子結合。相反，（1 1 1ˉ）晶面上的羥基銨陽離子顯露比較明顯，甲酸和水分子容易與（1 1 1ˉ）晶面上的羥銨陽離子相互作用。溶劑與晶面之間的相互作用能可以反映溶劑在晶面上的吸附能力。溶劑與晶面之間的吸附作用越強，晶面的相對生長速率越慢。相互作用能的計算方法見公式（5）。表5 為TKX?50 在5 種不同體積比的甲酸/水混合溶液中晶體習性的模擬結果。從表5 可以看到，溶劑與晶面之間的相互作用能都是負值，說明溶劑的吸附是放熱過程。溶劑與晶面之間的相互作用能與其結合能互為相反數。結合能越大，溶劑在晶面上的吸附越強。因此，相互作用能的絕對值越大，溶劑的吸附效果越強，使得晶面生長速率越慢。對于5 種體系而言，與混合溶劑產生較大相互作用能的晶面都是（1 1 1ˉ），產生較小相互作用能的晶面都是（0 2 0）。由此可以推測，不同體積比的甲酸/水混合溶劑都會阻礙（1 1 1ˉ）晶面的生長。

圖5 MD 模擬（0 1 1）和（1 1 1ˉ）晶面上甲酸/水的分子分布Fig.5 Formic acid/water molecular distributions of（0 1 1）and（1 1 1ˉ）crystal faces after MD simulation

前面根據公式（3）計算的S值推測，（1 1 0）晶面會與溶劑分子產生較強的吸附作用。S代表晶面的固有屬性，但是溶劑與晶面之間的相互作用的大小不僅取決于晶面的固有屬性，還取決于溶劑分子結構、溶劑分子電荷以及晶面?溶劑之間的結合方式與強度。因此，盡管（1 1 0）晶面的S值最大，但（1 1 0）晶面與混合溶劑之間的相互作用并不一定最大。此外，（1 1 0）晶面均具有最小的附著能絕對值，因此，（1 1 0）晶面生長速率相對最慢。表5 中同時顯示模擬結果的晶面面積所占總面積的比例，當甲酸/水的體積比為1/4、1/3 和1/2 時，（0 2 0）、（1 0 0）和（0 1 1）晶面完全消失；當甲酸/水的體積比為1/1 和2/1 時，（0 2 0）和（0 1 1）晶面完全消失。對于5 種體系而言，（1 1 0）面均是面積最大的晶面，具有最強的形態(tài)重要性。前面分析（1 1 0）晶面可能與溶劑分子產生較強的相互作用，而強相互作用導致（1 1 0）晶面生長較慢，這里的模擬結果與前面的分析結果是一致的。

圖6 為甲酸/水混合溶劑中TKX?50 晶體形貌的預測結果?？梢钥吹剑敿姿?水體積比為2/1 時，TKX?50 晶體形貌為片狀。當體積比為其它值時，TKX?50 晶 體 形 貌 為 棱 形。Li 等［15］通 過 實 驗 得 到 了TKX?50 在甲酸/水混合溶劑中體積比為1/1，溫度為298.15 K 時，重結晶得到的晶體形貌，如圖7 所示。對比圖6d 和圖7 可知，通過理論預測的形貌與實驗得到的形貌基本一致。當甲酸/水的體積比分別為1/4、1/3、1/2、1/1 和2/1 時，TKX?50 晶體的縱橫比分別為4.11、3.28、2.91、3.37 和6.71。Chen 等［28］預 測 的TKX?50 在水、乙二醇、二甲基亞砜、乙醇和甲苯溶劑中的晶體縱橫比分別為3.33、3.05、3.12、3.13 和3.20。劉英哲等［14］預測的TKX?50 在甲醇、四氫呋喃、乙酸乙酯和三氯甲烷溶劑中的晶體縱橫比分別為2.97、4.20、3.74 和4.26。綜合比較來看，當溶劑為甲酸/水，體積比為1/2 時，所預測的TKX?50 晶體形貌相對更接近球形，晶體形貌相對較好。

圖6 TKX?50 在不同體積比甲酸/水混合溶劑中預測的晶體形貌Fig.6 Crystal morphology prediction of TKX?50 in formic acid/water mixed solvent with different volume ratios

圖7 TKX?50 在 甲 酸/水 混 合 溶 劑 體 積 比 為1/1，溫 度 為298.15K 時，重結晶得到的晶體形貌［15］Fig.7 Crystal morphology of TKX?50 obtained by recrystalli?zation from formic acid/water mixed solvent at temperature of 298.15K and volume ratio of 1/1

3.3 溫度對TKX?50 結晶形貌的影響

溫度是影響晶體生長的重要因素之一，升高溫度會使溶液中溶質分子的運動速率加快，更容易克服溶劑的脫附能壘，能促進晶體的生長。另外，溫度的變化也會引起溶液過飽和度的變化，從而間接影響晶體的生長習性。因此可以通過調節(jié)溫度來控制晶面的生長速率，進而影響晶體的形貌。選擇甲酸/水的體積比為1/2，溫度分別設置為298，308，318 K 和328 K，執(zhí)行MD 模擬，模擬 結果歸納在表6?？梢钥吹?，不論溫度如何變化，混合溶劑與晶面之間的相互作用能都為負值，這說明溶劑的吸附是放熱過程，且不隨著溫度變化而改變。對于這4 種不同溫度的體系來說，與溶劑產生較大相互作用能的晶面都是（1 1 1ˉ），產生較小相互作用能的晶面都是（0 2 0），這說明甲酸/水混合溶劑會阻礙（1 1 1ˉ）晶面的生長。另外，不論溫度如何變化，（1 1 0）晶面均具有最小的附著能絕對值，說明（1 1 0）晶面生長速率相對最慢。表6 同時顯示，對于這4 種不同溫度的體系來說，（0 2 0）、（1 0 0）和（0 1 1）晶面都完全消失，而（1 1 0）和（1 1 1ˉ）晶面成為重要的生長面。隨著溫度升高，（1 1 0）面的面積逐漸增大，而（1 1 1ˉ）面的面積逐漸減小。最終形貌預測結果見圖8，預測的TKX?50 晶體形貌均為棱形，隨著溫度從298 K 增大到328 K，縱橫比分別為2.91、3.28、3.72 和4.87。當溫度為298 K 時，預測的TKX?50 晶體形貌相對更接近球形。

表6 TKX?50 在不同溫度下的甲酸/水混合溶劑中的能量和習性Table 6 Energies and habits of TKX?50 in formic acid/water mixed solvents at different temperatures

圖8 TKX?50 在不同溫度下甲酸/水混合溶劑中預測的晶體形貌Fig.8 Crystal morphology predictions of TKX?50 in formic acid/water mixed solvent at different temperatures

3.4 徑向分布函數

徑向分布函數（RDF）定義為給定一個粒子的坐標，距離這個粒子為r時出現其他粒子的概率，反映了體系中粒子的聚集特性。一般來說，溶劑與晶面之間的分子相互作用包括氫鍵（＜3.1 ?），范德華力（3.1～5.0 ?）和靜電力（＞5.0 ?）［13，31］。當甲酸/水的體積比為1/2，溫度為298 K 時，選擇TKX?50 的（1 1 0）晶面作為研究對象，溶劑與晶面之間的RDF 分析結果曲線如圖9 所示。對于TKX?50 晶面上的O 原子和溶劑中的H 原子而言，RDF 的第一個峰出現在2.2 ? 的位置，說明（1 1 0）─O 和溶劑?H 之間存在著氫鍵；第二個尖銳的峰出現在6.6 ? 的位置，說明（1 1 0）─O 和溶劑─H 之間有靜電力存在，而沒有范德華力存在。對于TKX?50 晶面上的H 原子和溶劑中的O 原子而言，RDF分別在2.3，3.4，4.6 和6.6 ? 的位置出現較為尖銳的峰，說明（1 1 0）─H 和溶劑─O 之間同時存在氫鍵、范德華力和靜電力?？傮w而言，氫鍵、范德華力和靜電力同時存在于TKX?50 的（1 1 0）晶面與溶劑分子間。

圖9 （1 1 0）晶面與甲酸/水混合溶劑體系中氧原子和氫原子之間的徑向分布函數Fig.9 The RDFs between oxygen atoms and hydrogen atoms in（1 1 0）crystal face and formic acid/water mixed solvent system

4 結論

本研究利用AE 模型預測了TKX?50 的真空形貌，通過修正的AE 模型預測了TKX?50 晶體在甲酸/水中的晶體形貌，并且比較了甲酸/水體積比和溫度對TKX?50 晶體形貌的影響，得到結論如下：

（1）TKX?50 真 空 中 的 重 要 生 長 面 為（0 2 0）、（1 0 0）、（0 1 1）、（1 1 0）和（1 1 1ˉ），縱橫比為1.98。（0 1 1）面具有最強的形態(tài)重要性，（1 1 1ˉ）面的生長速率最快，趨向于消失。此外，分析了晶面的特性，結果表明，（1 1 0）晶面最粗糙，因此（1 1 0）晶面可能會與溶劑分子產生較強的吸附作用。

（2）甲酸/水混合溶劑分子對TKX?50 各晶面的吸附強度不同，導致TKX?50 各晶面生長速率不同。當甲酸/水的體積比和溫度變化時，晶體形貌存在較為明顯的差異，然而其中共同點都是（1 1 0）晶面均具有最大的形態(tài)重要性。當甲酸/水的體積比為1/2，溫度為298 K 時，TKX?50 的晶體形貌相對更接近球形，縱橫比為2.91。

（3）徑向分布函數分析表明，TKX?50 的（1 1 0）晶面與溶劑分子間同時存在氫鍵、范德華力和靜電力。