N-丁基吡啶四氟硼酸盐/水二元体系的分子动力学模拟研究

王 玉,朱光来,王晨晨,徐建强,马赵鹏

(1.安徽师范大学 物理系 原子与分子物理研究所,芜湖 241002; 2.安徽达尔智能控制系统股份有限公司 研发中心,芜湖 241003)

1 引 言

室温离子液体是指在室温或接近室温下呈液态、由有机阳离子和无机阴离子所组成的离子化合物，简称为离子液体[1].离子液体具有良好的热稳定性和化学稳定性、极低的蒸汽压、可循环利用等优点[1-3].由不对称的阴、阳离子组成的结构特点，决定了离子液体具有这些独特的性质，研究离子液体的微观结构有助于了解离子液体的物理化学性质、溶剂特性等[3].

实验与理论研究告诉我们在离子液体及其溶液内部呈现静态或动态的非均质结构，存在着不同形态聚集体[4-7].离子液体与分子溶剂的混合体系中，组分比例不同会引起溶剂的热力学性质和局部微观结构的变化，对溶质分子的扩散、迁移等输运性质产生重要影响.而水作为生活中的常见分子溶剂，研究它与离子液体混合物的性质和结构具有重要的科学意义和应用价值[8].分子模拟方面，目前报道较多的是对烷基咪唑型离子液体与乙腈、水等混合物的径向分布函数、扩散系数、黏度的研究[9，10].对于吡啶离子液体，其热力学性质方面的实验研究已有较多报道[11-13]，而其与分子溶剂的混合物的性质与结构的分子模拟研究还有待进一步扩展[14，15].

本文以N-丁基吡啶四氟硼酸盐[BPy]BF4与水作为研究对象，用分子动力学模拟方法研究不同比例下二元体系的密度和不同组分的扩散系数，通过分析各组分间的径向分布函数、空间分布函数以及氢键作用，重点研究二元体系的微观结构.

2 力场与模拟方法

2.1 力场参数

[BPy]BF4和水分子结构如图1示意，其力场采用Lopes等人在OPLS_AA/AMBER力场基础上发展起来的全原子力场模型[16]，函数形式如下：

(1)

图1 [BPy]BF4和水分子结构示意图

其中Vtotal表示系统的总能量，包含了键长、键角、二面角、van der Waals(VDW)作用以及静电作用五个部分.其中非键作用力由Lennard-Jones (L-J)项和库仑项组成.键拉伸和键弯曲以简谐势表示.水分子的力场参数采用TIP4P四点结构模型[17].

2.2 模拟方法

分子动力学模拟采用GROMACS软件包[18]，结果分析也基于GROMACS提供的分析手段.模拟过程采用NPT系综，温度和压力维持在298 K和0.1 MPa.模拟步长为2 fs，每0.2 ps取样一次，范德华相互作用的截断半径以及库仑相互作用的截断半径均取为1.5 nm.模拟过程采用标准的周期性边界条件，用V-rescale方法维持体系温度，用Parrinello-Rahman方法维持体系压力，远程静电作用采用PME方法，键长约束采用LINCS算法.每个体系均先计算2 ns以使得体系构象达到平衡，然后在此基础上再运行4 ns用来进行统计分析.

3 结果与讨论

3.1 密度与自扩散系数

如表1 所示，二元体系的密度随着离子液体比例的增加呈递增的规律，模拟的数值与实验值(纯水0.997 g/cm3，纯离子液体1.2141 g/cm3)相比较[11]，误差均在允许的范围内.

表1 模拟得到体系的密度以及各组分的扩散系数

体系的自扩散系数D由Einstein关系(式(2))拟合得到：

(2)

式中t为时间，ri为组分i质心的矢量位移，本工作中各组分的自扩散系数D是通过线性拟合系统平衡后1ns时间内的均方位移(MSD)得到.从模拟结果中不难看出，随着离子液体比例的增大系统各组分的扩散系数均有明显的递减变化，这说明[BPy]BF4对该混合体系的扩散性质有很大的影响.另一方面来看，混合体系各组分的扩散系数随着水分子比例的减少而减小，这一现象预示着水分子的存在可能会有效降低[BPy]BF4离子液体的黏度.

3.2 微观结构

3.2.1径向分布函数

溶液的局部组成和微观分布情况可以通过分析它们的径向分布函数(RDFs)来描述，该函数的物理意义是指以某个粒子为中心，在距离为r处另一粒子的分布几率.图2～图6分别给出了二元体系中不同组分间的径向分布函数随离子液体摩尔分数xi的变化情况.

图2 阳离子上丁基链末端碳原子CT4之间的径向分布函数

图2表示阳离子丁基链末端碳原子(CT4)间的径向分布函数，其峰值随着离子液体比例的增加并没有明显的变化；图3表示阳离子吡啶环的HA1与阴离子的F原子间的径向分布函数，其峰值呈现出随离子液体比例的增加递增的规律.上述结果表明离子液体的浓度对阳离子丁基链组成的非极性区域的聚集程度没有多大影响，但对于吡啶环与阴离子组成的极性区域的影响却很大.这主要由于极性区域的主要作用是阴、阳离子间较强的库仑作用，而非极性区域则是烷基链之间的VDWs相互作用.

图3 阳离子吡啶环上HA1原子与阴离子F原子间的径向分布函数

图4和图5分别表示水分子中HW与HA1、F原子间的径向分布函数，虽然二者均随着离子液体比例的增加而递增，但图4中HW-F的峰值更高，且变化趋势明显，这表明混合体系中水分子与阴离子间的相互作用更强些，水分子更易靠近阴离子，原因可能是由于水分子与阴离子间的由氢键作用形成了网状聚集体.

图4 水分子HW原子与阴离子F原子间的径向分布函数

图5 水分子HW原子与阳离子吡啶环HA1原子间的径向分布函数

图6所示的是水分子中OW原子间的径向分布函数，很高的峰值说明水分子因较强的氢键作用聚集.而且其峰值随离子液体比例的增加也呈递增的趋势，可能是由于离子液体增加导致体系的黏度增大，扩散变慢，促使由氢键作用而形成的水分子小团簇更加稳定.

图6 水分子OW原子间的径向分布函数

3.2.2空间分布函数

空间分布函数(SDFs)能够更加直观、清晰地反映混合体系中各组分的空间分布情况.图7和图8分别表示离子液体的摩尔分数为0.1和0.5时的二元体系中不同组分间的空间分布函数.

图7 离子液体摩尔分数为0.1的混合体系中不同组分间的空间分布函数

图8 离子液体摩尔分数为0.5的混合体系中不同组分间的空间分布函数

由图7(a)和图8(a)可以看到，离子液体摩尔分数xi为0.1和0.5时，阴离子周围的水分子分布几率图形非常相似，在四个F原子周围的都近似呈均等分布，这是由于F原子与水分子之间因OW-HW…F氢键作用形成了聚集体，但xi为0.5时的isovalue值要高许多，说明这时水分子在阴离子周围的几率更大，这一结果与HW-F原子间的径向分布函数随离子液体比例增大而增加的规律相一致.由图7(b)和图8(b)可以看到，在xi为0.1时，阳离子周围的水是环绕整个阳离子头部分布的，这是因为水分子可以与阳离子吡啶环上的氢原子形成不同强度的氢键，而xi为0.5时，水主要分布在一侧，此时水在阴离子周围的分布更有序；如图7(c)和图8(c)所示，阴离子在阳离子周围主要分布在头部，这主要是由于阴、阳离子之间存在较强的库仑作用所致.

3.2.3氢键作用

通过观察图3和图4所示的F-HA1、F-HW之间的径向分布函数可以发现，两者第一个波峰的位置在0.170-0.253 nm范围内，均短于H、F原子间的Van der Waals半径之和(0.267 nm)[19]，这表明阳离子与阴离子、水分子与阴离子之间分别有类似C-HA1…F和OW-HW…F氢键形成.进一步比较它们的峰位可以看出，OW-HW…F氢键(最大峰对应位置约为0.172 nm)要强于C-HA1…F氢键(约为0.253 nm)，并且这两种氢键作用都随着水分子含量的减少而增强.

为了探究水分子自身之间形成的氢键性质随离子液体变化的趋势，我们计算了不同离子液体比例下二元体系中每个水分子周围形成的平均氢键数目以及氢键的存活寿命，如图9所示.由图可知，随离子液体的增加，水分子数目明显减少，体系中每个水分子周围形成的平均氢键数目减少，而氢键寿命则增大，这与文献报道的水分子氢键的生存周期随溶质浓度增大而增大的规律类似[20].此时离子液体各组分的扩散系数呈递减的趋势，说明离子液体比例越高各组分扩散的速度越慢，使得水分子间的运动也变得缓慢，增加了氢键断裂的难度，其氢键寿命得以延长.

图9 水分子的平均氢键数目与氢键寿命随离子液体比例的变化

4 结 论

利用分子动力学模拟方法研究分析了吡啶离子液体[BPy]BF4与水组成的二元体系的微观结构和各组分之间的相互作用，主要结论如下：水的存在会有效地降低离子液体的黏度，各组分的自扩散系数随着水分子比例的减少而变小；随着离子液体比例的增大，除了丁基链末端碳原间的径向分布函数无明显变化外，其它组分原子间的径向分布函数的峰值均随之增加 说明水与阴离子、水与阳离子吡啶环之间的相互作用逐渐增强，并且水与阴离子间的相互作用更强一些，水分子在微观分布上更容易趋近阴离子，这与三者之间形成的氢键有密切关系；空间分布函数则直观地反映出阴离子周围的水的分布几率近似呈均等分布，阴离子主要分布在阳离子的头部周围，而且水分子在阴离子周围的分布几率随离子液体比例增大而增加；通过分析水分子间的氢键还发现，每个水分子周围的平均氢键数目随离子液体比例的增加而减少，而氢键寿命则随之增大.