TKX⁃50 在甲酸/水混合溶剂中生长形貌的分子动力学模拟

周 涛，陈 芳，李 军，曹端林，王建龙

（中北大学化学工程与技术学院，山西 太原 030051）

1 引言

5，5′‑联四唑‑1，1′‑二氧二羟铵（TKX‑50）是一种四唑类含能离子盐，区别于传统的含能材料，TKX‑50分子结构中并没有硝基致爆基团，但是爆速和爆压等爆轰性能均超过了奥克托今（HMX）［1］，主要原因是环形骨架中存在高能的N—N、N=N 键。自从TKX‑50被合成后，国内外针对其衍生物设计以及合成工艺路线优化进行了广泛的研究［2-6］。一般来说，含能材料的晶体形貌是影响其性能的重要因素之一［7］，在实际应用中经常期望得到特定的形貌，以满足其安全性能需求。例如，晶体形貌规则且接近球形，可以提高炸药的装药密度和流散性，进而改善药柱的力学性能和机械感度。徐容等［8］利用溶剂‑非溶剂重结晶法得到粒度跨 度 小、撞 击 感 度 低 的2，6‑二 氨 基‑3，5‑二 硝 基 吡嗪‑1‑氧化物（LLM‑105）晶体；周诚等［9］在二甲基亚砜/水（DMSO/H2O）混合溶剂中重结晶得到表面光滑、热感度较低的1，1‑二氨基‑2，2‑二硝基乙烯（FOX‑7）晶体；许诚等［10］采用降温重结晶法获得粒度跨度小、表面光滑的TKX‑50 晶体，降低了其机械感度。

近年来，计算机科学的蓬勃发展，使得人们逐渐掌握了在原子分子水平上模拟晶体生长过程的手段。其中，分子动力学（MD）方法已经广泛地被应用于研究炸药晶体生长过程中的溶剂效应。Duan 等［11］借助MD 方法解释了丙酮溶剂对HMX 晶体形貌的影响。Chen 等［12］通过MD 方法计算了黑索今（RDX）晶体在丙酮溶剂中的生长形貌；Lan 等［13］采用MD 方法模拟了外部生长环境对六硝基六氮杂异伍兹烷（HNIW）晶体形貌的影响，外部生长环境包括溶剂、温度和过饱和度等。关于TKX‑50 在一元溶剂中的晶形预测已有报道［14］，但是TKX‑50 在二元混合溶剂中的晶体形貌尚未被系统研究，混合溶剂体积比与温度对TKX‑50 晶体形貌的影响情况尚不清楚。而TKX‑50 在大多数常用溶剂中的溶解度非常差，如二甲基甲酰胺（DMF），甲醇，乙醇，甲苯，乙酸乙酯等，但在甲酸和水中具有更好的溶解性［15］。因此，本研究采用MD 方法模拟TKX‑50在甲酸/水混合溶剂中的晶体形貌，采用附着能（AE）模型预测TKX‑50 在真空中的晶体形貌，预测TKX‑50 在不同温度、不同体积比的甲酸/水中的结晶形貌，并通过实验结果来验证模拟的准确性；最后进行径向分布函数分析，以期了解TKX‑50 晶面与混合溶剂之间存在的相互作用力。

2 计算细节

2.1 AE 模型及其修正

AE 模型是在周期性键链（PBC）理论基础上建立起来的模型［16-17］，它根据晶体对称性和分子间键链性质计算晶体的附着能。该理论认为各个晶面相对生长速率与晶面附着能的绝对值成正比［18］：

式中，Rhkl为相对生长速率；Eatt为晶面的附着能，kJ·mol-1。

晶面附着能绝对值越高，晶面生长速率越快，晶面生长趋于减小或消失；晶面附着能绝对值越低，晶面生长速率越慢，在最终的晶体形态中越容易得到显露。通过评估晶体不同晶面的相对生长速率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 的晶体形貌，确定形态学上重要的晶面。然后对TKX‑50 的重要晶面进行切割，将其构造为3×3×3 的周期性超晶胞结构模型，并对超晶胞结构模型进行优化。随后，构建溶剂层模型，溶剂体系为不同体积比的甲酸/水（1/4，1/3，1/2，1/1 和2/1）混合溶剂。对溶剂层进行几何优化，将优化后的溶剂层沿着c轴与晶面层对接，构建晶体‑溶剂双层模型，以此来探究溶剂对晶体形貌的影响。晶面层与溶剂层之间的真空距离为3 Å，而溶剂层上方空出50 Å 的真空距离，以便消除上下晶体表面的影响。对双层模型进行20000 次的迭代优化，然后进行MD 模拟，MD 模拟期间要固定晶面层。MD 模拟总时间设定为300.0 ps（300000 fs），时间步长为1 fs，系综选择NVT，温度由Andersen 恒温器控制［25］。当温度和能量的波动小于10%时，整个体系可以认为达到了平衡。

3 结果与讨论

3.1 真空中TKX⁃50 的结晶形貌

TKX‑50 晶体单晶衍射数据与优化后的晶胞参数列于表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 晶胞参数的实验数据和优化数据的对比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 真空中晶体晶形为不规则长块形，主要包括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ˉ），晶体晶形为不规则长块形。本研究模拟的结果与文献［14］中的模拟结果相一致。纵横比定义为晶体习性的最长和最短直径之间的比值。纵横比越接近1，预测的晶体形貌越接近球形。本研究模拟预测TKX‑50 真空形貌的纵横比为1.98，其形貌参数列于表3。其中，（0 1 1）面的面积比最大，达到了35.65%，具有最大的形态重要性。（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 形态重要晶面的参数Table 3 The parameters of morphologically important faces of TKX‑50 in vacuum

TKX‑50 分子的羟铵阳离子和联四唑阴离子通过强静电作用结合，不同的分子排布方式可能会导致不同的电荷分布，最终影响晶面与溶剂之间的相互作用。TKX‑50 各个重要晶面的分子堆积结构见图4。可以看到，在（0 2 0）晶面上有排列规则的联四唑阴离子显露。（1 0 0）晶面上联四唑阴离子环对立分布，空间位阻较大。（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）晶面由于其表面的结构特征，可以对溶质分子和溶剂分子产生较强的吸附能力。但是，晶面更容易吸附溶质分子还是溶剂分子，仅分析表面结构是不够的。溶液结晶是一个复杂的相转变过程。一方面，溶剂吸附于晶体降低了界面能，使晶面由光滑面转变为粗糙面，这会促进晶面的生长；另一方面，溶剂的优先吸附占据了晶面的生长活性位，溶质生长必须先克服溶剂的脱附能垒，这会阻碍晶面的生长［27］。溶质分子和溶剂分子在界面上竞争吸附，其吸附强弱能力需要从多方面来分析，比如晶面结构，溶质和溶剂分子结构，表面静电势等。本课题组前期工作 分 析 了TKX‑50 晶 面 的 表 面 静 电 势［28］，认 为 除 了（1 0 0）晶面，其它晶面都为极性面。晶面的正电荷比较多，因此溶质中的阴离子相较于阳离子来说更容易吸附到晶面上。此外，选择具有负电子基团或富电子芳环的溶剂分子，更容易控制TKX‑50 晶面的生长速率，从而控制TKX‑50 的晶体形貌。

3.2 不同体积比的甲酸/水混合溶剂对TKX⁃50 晶体形貌的影响

溶质分子从溶液主体扩散到晶体表面并且长入晶体，促进了晶体的生长。然而，溶剂会优先占据晶面的生长活性位点，溶质分子吸附到晶面必须克服溶剂的脱附能垒［29-30］。因此，溶剂在晶面上的吸附会影响晶面的生长速率，促使不同晶面的各向异性生长，导致晶体的生长形貌不同。以（0 1 1）和（1 1 1ˉ）晶面为例，MD 模拟（0 1 1）和（1 1 1ˉ）晶面上甲酸/水的分子分布如图5 所示。可以看到，（0 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）面均是面积最大的晶面，具有最强的形态重要性。前面分析（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 结晶形貌的影响

温度是影响晶体生长的重要因素之一，升高温度会使溶液中溶质分子的运动速率加快，更容易克服溶剂的脱附能垒，能促进晶体的生长。另外，温度的变化也会引起溶液过饱和度的变化，从而间接影响晶体的生长习性。因此可以通过调节温度来控制晶面的生长速率，进而影响晶体的形貌。选择甲酸/水的体积比为1/2，温度分别设置为298，308，318 K 和328 K，执行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 之间同时存在氢键、范德华力和静电力。总体而言，氢键、范德华力和静电力同时存在于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）面具有最强的形态重要性，（1 1 1ˉ）面的生长速率最快，趋向于消失。此外，分析了晶面的特性，结果表明，（1 1 0）晶面最粗糙，因此（1 1 0）晶面可能会与溶剂分子产生较强的吸附作用。

（2）甲酸/水混合溶剂分子对TKX‑50 各晶面的吸附强度不同，导致TKX‑50 各晶面生长速率不同。当甲酸/水的体积比和温度变化时，晶体形貌存在较为明显的差异，然而其中共同点都是（1 1 0）晶面均具有最大的形态重要性。当甲酸/水的体积比为1/2，温度为298 K 时，TKX‑50 的晶体形貌相对更接近球形，纵横比为2.91。

（3）径向分布函数分析表明，TKX‑50 的（1 1 0）晶面与溶剂分子间同时存在氢键、范德华力和静电力。