劉英哲

(西安近代化學研究所，陜西 西安 710065)

引 言

晶形控制是目前含能材料改性研究的重要途徑之一[1]。對于相同的含能晶體，晶形能夠在很大程度上影響感度、能量、流散性等性質[2-3]。例如，在相同尺寸下，球形晶體通常比其他晶形具有更高的裝填密度和更低的機械感度[4-5]。然而，含能材料的晶形控制仍然較為困難，因為晶形生長不僅與晶體結構有關，還受到溶劑、溫度、添加劑、結晶工藝等多種條件的影響。溶液結晶是含能材料最常用的晶形生長方式，但采用不同的溶劑時，得到的晶形差異很大[6]。

與基于“試錯法”的結晶實驗相比，通過理論方法預測晶形，有助于從微觀上深入理解晶體生長機制，能夠快速篩選溶劑、溫度等結晶條件，從而為含能材料的晶形控制提供理論指導。當前，含能材料的晶形預測方法主要包括：Bravais-Friedel-Donnay-Harker (BFDH) 準則[7]、附著能模型[8-9]、占有率模型[10-11]、蒙特卡洛模擬方法[12-14]、界面結構分析模型[15]、螺旋生長模型和二維成核模型[16-19]等，這些方法的原理逐漸從結構因素、能量因素向生長機制演化，預測精度也逐漸增加。其中，附著能模型在考慮晶體結構因素的基礎上，又引入了能量因素，預測結果相對準確，計算簡單，又有很好的軟件支持，成為含能材料領域中應用最廣泛的晶形預測方法。

為此，本文綜述了國內外近十年來附著能模型在含能材料中的應用情況，詳細闡述了附著能的理論基礎與計算方法，討論了影響晶形預測結果的重要因素，并介紹了附著能模型的最新研究進展。

1 附著能模型理論與方法

1.1 基礎理論

附著能模型基于恒定生長速率的簡單假定[8-9]，即晶體中不同晶面的相對生長速率(Rhkl)與其附著能(Eatt)的絕對值成正比關系：

Rhkl∝|Eatt|

(1)

晶面的附著能絕對值越大，其相對生長速率越快，在生長過程中越容易消失。附著能定義為一層厚度為dhkl的晶片附著在晶面(h k l)上所釋放的能量，可由式(2)進行計算：

Eatt=Elatt-Eslice

(2)

式中：Elatt為晶格能量；Eslice為晶片能量。對于穩(wěn)定生長的晶面，附著能為負值。需要注意的是，不同晶面的Eatt需要在相同的化學計量下(如每分子、每晶胞)比較才有意義，目前在絕大部分文獻中，Eatt的單位簡寫為kcal/mol，實際上完整的單位應為kcal/(mol·molecule)或kcal/(mol·unitcell)。

基于給定的晶體結構，通過計算不同生長晶面的附著能，可預測真空條件下的生長形貌。在溶劑條件中，溶劑將附著在生長晶面上，溶質在晶面上生長時要克服溶劑層的阻力，導致晶面生長速率受到抑制。為此，需要對附著能進行修正：

(3)

式中：Es為溶劑吸附所引入的能量修正項，表示溶劑對晶面生長的抑制作用；S為修正因子，體現(xiàn)了生長晶面的結構特征，可理解為不同晶面結構對Es的權重，其定義為:

(4)

式中：Aacc為晶面單元(h k l)的溶劑可及表面積；Ahkl為晶面單元的橫截面積。顯然，S值越大，晶面越粗糙；S值越小，晶面越平坦。Es可通過溶劑與晶面的相互作用能(Eint)進行計算：

(5)

式中：Abox為模擬盒子中生長晶面(h k l)的橫截面積；Eint為晶面-溶劑體系的總能量(Etot)與晶面能量(Ecry)和溶劑能量(Esol)之差，即：

Eint=Etot-(Ecry+Esol)

(6)

在溶劑條件下，假定晶面的相對生長速率與修正后的附著能仍然成正比關系，即：

(7)

因此，通過計算溶劑中的修正附著能，即可預測晶體在溶劑中的生長形貌。為了便于理解溶劑條件下附著能公式中的相關計算參數(shù)，將其物理含義形象地表示出來，如圖1所示。另外，基于類似的思想，也可將附著能模型進一步拓展至添加劑條件下[20-22]，此處不再贅述。

圖1 溶劑中附著能計算公式示意圖Fig.1 Schematic diagram of attachment energy calculation in solvent

1.2 計算模型

為了預測溶劑中的生長形貌，需要構建合理的晶面-溶劑界面模型，從而準確地計算附著能。圖2為晶面-溶劑界面模型的構建示意圖。首先，基于給定的晶胞結構，計算低米勒指數(shù)晶面的附著能，得到真空中的預測晶形和穩(wěn)定生長晶面。然后，沿著每個生長晶面的(hkl)指數(shù)，對晶胞進行切割獲得生長晶面單元，并構建超晶胞，從而獲得三維的晶面模型。對于溶劑模型，可根據(jù)給定密度將溶劑分子隨機地填充到三維盒子中。最后，沿著生長晶面的法向量方向，將溶劑模型“拼接”到晶面模型上方，形成晶面-溶劑界面的初始結構。在添加劑條件下，也可通過類似方法構建計算模型，但一般用單個添加劑分子替代溶劑層結構[23]。

圖2 晶面-溶劑界面模型搭建示意圖 Fig.2 Scheme in construction of crystal-solvent interfacial model

隨后，在NVT系綜中對上述界面結構進行分子動力學模擬。由于施加了三維周期性邊界條件，為了盡可能地避免溶劑與晶面的鏡像發(fā)生相互作用，還需要在溶劑層上方設置一定厚度的真空層。真空層的厚度要大于非鍵作用的截斷半徑，通常為幾十埃。另外，附著能模型假定生長晶面上不發(fā)生結構重排[8]，一般可在分子動力學模擬中固定晶面的位置。待界面體系到達平衡后，對模擬軌跡進行平均求得晶面-溶劑相互作用能，從而進一步計算溶劑中的附著能。

2 附著能模型在含能材料中的應用

2.1 在晶形預測中的應用

當前，附著能模型已經廣泛應用于含能材料的晶形預測，研究體系包括常見的分子晶體、離子晶體和含能共晶，分別見表1、表2和表3。這些研究預測了含能材料在真空、(混合)溶劑、添加劑等不同生長條件的晶形，并討論了溫度效應對晶形生長的影響。此外，將附著能結果進一步結合徑向分布函數(shù)、擴散系數(shù)等計算分析方法，從微觀上揭示了溶劑與生長晶面的相互作用，為含能材料的晶形控制，甚至是共晶形成提供了理論支撐[24-26]。

表1 含能分子晶體的晶形預測文獻匯總Table 1 Summary on the morphology predictions of energetic molecular crystals

表2 含能離子晶體的晶形預測文獻匯總Table 2 Summary on the morphology predictions of energetic ionic crystals

表3 含能共晶的晶形預測文獻匯總Table 3 Summary on the morphology predictions of energetic cocrystals

2.2 不同計算結果間的差異

對于大部分含能晶體，附著能模型能夠從定性上給出與實驗數(shù)據(jù)相吻合的預測結果。然而，從相同體系的預測晶形預測對比圖(見圖3)可以看出，有時附著模型預測的晶形也存在一定差異。盡管附著能模型的理論假設非常明晰，但附著能的計算結果受很多因素影響，如力場參數(shù)，模型尺寸，晶面-溶劑相互作用能等，這都將直接決定晶形的預測結果。因此，如何準確地計算附著能是一個關鍵問題。

圖3 相同體系下的晶形預測結果對比Fig.3 Comparison on the morphology predictions of same systems

3 附著能的準確計算

3.1 Eatt的計算

從附著能的定義上看，附著能的本質是相互作用能。因此，在實際計算過程中，真空中的附著能可通過式(2)求得：

(8)

為了準確計算附著能，需要考慮其隨著晶面厚度的收斂情況。如圖4所示，可計算一系列不同厚度的晶面能量，當其能量差趨于穩(wěn)定時，則附著能達到收斂。此處以β-HMX晶體為例，晶胞的Z值為2，5個生長晶面(0 1 1)、(1 1 0)、(0 2 0)、(1 0 1)和(1 0 -1)的Z值分別為2、2、1、2和2。

圖4 真空中附著能計算示意圖 Fig.4 Schematic diagram of attachment energy calculation in vacuum

表4給出了真空中HMX生長晶面的附著能隨晶面厚度的變化情況[31]。不同晶面的附著能對晶面厚度的收斂情況是不同的。顯然，當晶層厚度不足時，附著能計算偏差較大。對于HMX晶體，當n≥6時，所有晶面的附著能都已收斂，與Materials Studio軟件[68]給出的結果基本一致。另外，考察附著能隨晶面厚度的收斂情況，也可為界面模型中晶面層厚度的構建提供參考。

表4 真空中HMX生長晶面的附著能隨晶面厚度的變化情況Table 4 Attachment energies of HMX growth faces in a vacuum as a function of crystal layer thickness (unit: kcal/mol/unitcell)

3.2 力場參數(shù)

附著能的計算是基于分子力場完成的。因此，在開始計算附著能之前，要對力場參數(shù)的可靠性進行評價。選擇合適的分子力場是準確計算附著能的基礎。通常，可將擬選用的分子力場用于晶體結構優(yōu)化，根據(jù)優(yōu)化后的晶格參數(shù)與實驗數(shù)據(jù)的偏差，對力場的準確性進行驗證。此外，也可基于分子力場計算晶格能，與其實驗值進行對比。晶格能可通過實驗測定的升華焓(ΔHsub)數(shù)據(jù)進行求解[8,29]：

Elatt=-ΔHsub-2RT

(9)

式中：R和T分別為氣體常數(shù)和溫度；2RT為校正因子。晶格參數(shù)與晶格能考察的主要是分子間作用參數(shù)，同時也要注意晶體優(yōu)化后分子結構是否發(fā)生較大變化。另外，可通過溶劑密度計算來驗證分子力場對溶劑的兼容性。

COMPASS力場[69]是目前含能材料領域中應用最廣泛的力場，可適用于絕大部分CHON含能分子。然而，隨著近年來新型含能化合物的陸續(xù)合成，有一些特殊的成鍵類型，難以用COMPASS力場描述，比如TKX-50中的N→O鍵。表5對比了COMPASS、PCFF[70-72]、CVFF[73]3種力場對TKX-50晶格參數(shù)的優(yōu)化結果[23]。從表5可知，COMPASS力場優(yōu)化的晶格參數(shù)與實驗數(shù)據(jù)相差較大。CVFF力場優(yōu)化的b軸長度與實驗值也存在較大偏差。相比之下，只有PCFF力場的優(yōu)化結果與實驗結果較為接近，因此被用于描述TKX-50晶體[21,59,60]。

表5 不同力場優(yōu)化下的TKX-50晶格參數(shù)Table 5 Lattice parameters of TKX-50 crystal optimized by different force fields

需要注意的是，TKX-50是一種含能離子鹽，正負離子之間存在很強的靜電相互作用，對生長晶面的附著能影響很大。在PCFF力場中，TKX-50的正、負離子分別帶有+1e和-2e的電量，與量化計算結果+0.82e和-1.64e不相符[74]。因此，可通過密度泛函理論GGA-PBE/DNP對TKX-50晶胞進行電荷布居分析，將Mulliken電荷作為原子電荷，進一步修正TKX-50晶體內的靜電作用。經該參數(shù)優(yōu)化后，TKX-50的晶格參數(shù)變?yōu)閍=5.45?,b=10.89?,c=6.23?,α=90.00°,β= 8.44°,γ=90.00°，與實驗數(shù)據(jù)也較為相近。

表6匯總了真空下TKX-50的晶形預測結果。同樣采用PCFF力場描述TKX-50時，晶形預測結果也略有不同，但主要生長晶面和長徑比基本相同。對PCFF力場的電荷參數(shù)修正后，晶形與生長晶面發(fā)生了較大變化。另外，COMPASS力場計算的晶形類似于紡錘狀。由此可見，晶形預測結果對力場參數(shù)較為敏感。

除了TKX-50，劉寧等[44]對FOX-7的Dreiding力場參數(shù)[75]進行了部分修正，使晶格參數(shù)的相對計算誤差在5%之內；Song等[40]結合COMPASS力場和Gasteiger電荷，進一步將晶格參數(shù)的相對計算誤差控制在3%以內，均比COMPASS力場(7%)有了一定提升[41]。另外，Song等[76]結合COMPASS力場參數(shù)和RESP電荷描述DNTF晶體，晶胞參數(shù)的最大相對誤差僅為-2.86%，相比之下，采用COMPASS力場中的默認電荷時，晶胞參數(shù)的最大相對誤差達-6.00%[77]。

3.3 模型尺寸

在分子動力學模擬中，模型尺寸也是一個影響附著能計算結果的重要因素。顯然，模型越大，附著能計算越精確，耗時越長。相反，模型尺寸也不能過小，至少要大于非鍵作用截斷半徑的兩倍，以避免周期性邊界條件下分子與其鏡像發(fā)生相互作用[78]。因此，如何搭建合適的模型體系，在準確性與耗時性上取得平衡，是一個需要考慮的問題。在構建晶面-溶劑界面模型時，研究人員通常采用3×3×3超晶胞的晶面模型和幾百個分子的溶劑模型，并沒有一個通用的構建規(guī)則。

為此，Li lijie等[30,36,37,41,46]系統(tǒng)研究了模型尺寸對附著能計算結果的影響，包括晶面層、溶劑層、真空層的厚度以及模擬盒子的長度和寬度。首先，考察了非鍵作用能隨截斷半徑的收斂情況，認為非鍵作用的截斷半徑(dc)選取15.5?較為合適。然后，分別計算了溶液中附著能隨晶面層、溶劑層、真空層厚度和模擬盒子長度、寬度的變化。最后，通過對β-HMX、FOX-7、ε-CL-20、NTO多種含能材料的測試，給出了一個合適尺寸的模型構建規(guī)則：晶面層的長度、寬度和厚度應大于2dc、2dc和dc，溶劑層和真空層的厚度應大于dc，即模型尺寸應在圖5填色區(qū)域之外。基于該模型尺寸，溶液中生長形貌的預測結果與電鏡數(shù)據(jù)極為吻合。

圖5 模擬盒子的建議尺寸 Fig.5 Suggested size of simulation box

3.4 Eint的計算

探索溶劑在晶面上的微觀吸附行為，有助于理解溶劑效應對晶面生長的影響，從而更加準確地計算晶面-溶劑相互作用。研究表明[24,31,63,67,79-82]，在晶面-溶劑作用的影響下，溶劑在晶面上呈現(xiàn)出了典型的結構特征——層狀分布。此處，以HMX和丙酮溶劑為例，進一步闡述溶劑的吸附特點。

圖6展示了丙酮溶劑(紅色)和HMX(黑色)沿生長晶面(0 0 1)法向的質量密度分布情況[31]。顯然，溶劑分子在HMX晶面法向上并不是均勻分布，而是呈層狀分布。隨著溶劑分子與HMX晶面法向距離的增加，溶劑峰的波動逐漸減小，直至局域密度趨于體密度，說明此處溶劑分子已經不受晶面作用的影響，與體相中的分布一致。因此，若要完整體現(xiàn)出溶劑在生長晶面上的結構排布，溶劑層需要具有一定的厚度。盡管不同晶面結構會導致局域密度有所不同，但均呈現(xiàn)出層狀分布的特征。

圖6 沿晶面法向的質量密度分布Fig.6 Mass density profiles along the normal to the crystal face

為了分析層狀分布對晶面-溶劑相互作用能的影響，根據(jù)溶劑峰的個數(shù)將整個溶劑層進行分層，如圖6中為L1、L2和L3的色塊，然后分別計算區(qū)域L1、L1+L2、L1+L2+L3和整個溶劑層與生長晶面的相互作用能。對于HMX和丙酮體系，HMX晶面與溶劑層L1的作用能約占總相互作用能的90%以上?；诓煌軇拥腅int，最終預測的晶形結果較為相似，只是長徑比略有不同[31]。

需要注意的是，在動力學模擬過程中，個別溶劑分子由于自由擴散作用可能會吸附在生長晶面的鏡像上，為了確保計算精度，在計算Eint時要將這些溶劑分子剔除。如果采用了“三明治”計算模型，即生長晶面上下均有溶劑層，則不存在這個問題[76]。同時，還應注意附著能的化學計量。在真空中計算附著能時，若以晶胞為計量單位，需要將能量校正項Es轉換為與Eatt相同的單位。

4 附著能模型的發(fā)展

4.1 溶劑效應的校正

如何更加合理準確地描述溶劑效應對生長形貌的影響是附著能模型發(fā)展的一個方向。從式(2)可知，能量校正項Es體現(xiàn)了晶面與溶劑的相互作用，校正因子S相當于生長晶面結構帶來的權重。值得注意的是，Eatt和Es的本質都是相互作用能，區(qū)別在于：前者是晶面與生長晶片的相互作用，而后者反映了晶面與溶劑的相互作用。然而，在真空中計算Eatt時，并沒有用S進行校正。因此，在溶液中計算附著能時，也不應采用S對Es進行校正。這是因為：在計算相互作用能時，已經隱含了晶面結構的影響，即晶面越粗糙，相互作用的接觸面積也就越大，同樣條件下的相互作用能也就越強，無需進一步修正。

為了更好地解釋說明，定義了接觸面積：

(10)

其中，Atot、A1和A2分別為體系和所包含組分的溶劑可及表面積。以HMX為例，計算表明，在真空中和溶劑中，晶面分別與生長晶片和溶劑的接觸面積幾乎是一樣的[31]。因此，在同樣的接觸面積下，生長晶面在真空中的相互作用能Eatt和在溶劑中的相互作用能Es的計算標準應是一致的。實際上，校正因子S對晶面-溶劑相互作用額外加了權重，將對溶劑中的附著能帶來一定偏差，特別是結構相對粗糙的生長晶面，有時可能使得附著能為正值，致使晶面不穩(wěn)定。最近的研究基于這一新的校正公式，對溶劑中HMX的晶形進行了預測，計算結果與實驗結果吻合較好[30-31]。

4.2 過飽和度的探索

附著能模型假設晶面的生長速率與附著能成正比，即每個晶面的生長速率是恒定的。因此，從原理上講，附著能模型難以處理生長速率變化顯著的生長條件，如過飽和度等[7]。

最近，Li lijie等[30,37]在構建溶劑層時，將一定數(shù)量的溶質分子也添加在內，研究了過飽和度對β-HMX和ε-CL-20在溶劑中生長形貌的影響。過飽和度(δ)可通過式(11)在模型構建時進行計算：

(11)

式中：C為計算模型中溶質的實際濃度，由溶劑層中溶質和溶劑的數(shù)量決定；C0是實驗測定的飽和濃度。晶形預測結果表明，不同的過飽和度對形貌的長徑比影響較大。由于結晶實驗中的過飽和度很難保持不變，因此晶形預測結果尚不能與實驗結果進行對比。

4.3 晶形預測新策略

圖6表明溶劑沿生長晶面法向呈層狀分布，當溶劑分子遠離晶面一定距離后，其與晶面的作用可忽略不計，此時溶劑可看作自由擴散。當溶劑距離晶面較近時，溶劑吸附在晶面上，對溶質在晶面上的生長起到一定阻礙作用。

最近的分子動力學模擬表明，溶劑分子能夠長時間地停留在晶面上的某些特定位置(吸附位點)[82]，具有最強的結合能力，對溶質生長的阻礙作用最大，這說明溶劑與晶面的相互作用也是不均勻的。然而，在計算溶液中的附著能時，溶劑效應是通過“平均化”的晶面-溶劑相互作用能來描述的(見式(6))。這種溶劑層與晶面層的相互作用并不能真正體現(xiàn)出對溶質生長最為重要的抑制效應。

實際上，吸附位點處的結合能可以看作是“木桶效應”中的最短板，是決定晶體生長速率的關鍵因素。圖7給出了晶體表面生長過程的示意圖。由于溶劑在吸附位點處的結合能力最強，溶質將吸附位點處的溶劑替換則可以看作是晶面生長的“限速”步驟。為此，基于附著能模型的理論框架，提出了晶形預測新策略，即采用吸附位點處溶劑與溶質的結合能代替晶面-溶劑相互作用能，采用吸附位點處溶質與溶質的結合能代替真空中的附著能，具體計算細節(jié)可參考文獻[82]。顯然，相比于溶劑中原有附著能的計算，從晶體生長過程上看，這種策略更加合理?；谠摬呗?，預測了HMX在丙酮溶劑中的晶形，計算結果與實驗結果吻合良好[82]。

圖7 晶體表面生長過程示意圖 Fig.7 Scheme in the crystal growth process in solvent

5 結束語

附著能模型是目前含能材料領域中應用最為廣泛的晶形預測方法，能夠定性地給出溶劑、溫度等不同條件下的生長形貌，可為含能材料的晶形控制提供理論參考。在使用附著能模型時，需要充分考慮力場選擇、模型構建、溶劑效應等因素，從而準確地計算附著能。否則，相同體系的預測結果也可能會存在差異。值得注意的是，附著能模型僅僅基于恒定生長速率的簡單假設，尚不能考慮結晶實驗中的諸多因素，預測晶形有時與實驗結果偏差較大也是合理的。為此，仍需進一步發(fā)展基于生長機制的晶形預測方法。