聂文洁， 王剑飞， 赵贯甲， 尹建国， 马素霞

(1. 太原理工大学电气与动力工程学院， 太原 030024；2. 太原理工大学循环流化床高效清洁与利用山西省重点实验室， 太原 030024；3.太原锅炉集团有限公司， 太原 030024)

1 引 言

低共熔溶剂(deep eutectic solvents)是一种新型的绿色溶剂，不仅具备离子液体的蒸汽压低、溶解性好、导电性强、化学稳定性好等诸多优势，还克服了离子液体制备复杂、提纯困难、微毒、价格昂贵、不易降解等缺点[1]. 2003年，Abbott等[2]首次将氯化胆碱(choline chloride，ChCl)和尿素(urea)在1∶2的摩尔比下加热搅拌混合，冷却后形成了一种熔点很低的混合液体，即低共熔溶剂.由于低共熔溶剂制备简单、原料绿色廉价，且兼具离子液体良好的化学性质，引起了广大研究者的重视. 其在酸性气体吸收、分离过程、化学反应、生物催化和电化学等领域显示出良好的应用前景[3, 4].

低共熔溶剂至少由两种成分组成：一种称为氢键受体(HBA)，一种称为氢键供体(HBD).氢键受体通常有季铵盐(如氯化胆碱)和两性离子(如甜菜碱)等；氢键供体主要有尿素、醇类(如乙二醇、甘油、丙二醇、丁二醇)、酸类(乙酰丙酸、丙二酸、草酸)和苯酚等[1]. 由于每种HBA和HBD都能随意组合，因此理论上可形成无限多种低共熔溶剂. 其组成不同，结构和性质也不同. 氢键的形成是导致低共熔溶剂熔点降低和宏观热力学性质变化的重要因素. 分子动力学模拟提供了一种可以计算流体宏观特性和微观结构的方法，是研究这类体系经济、可靠的手段. Logotheti等人[5]采用分子动力学模拟计算了298.15～333.15 K温度下甲基咪唑类离子液体的密度、扩散系数和热膨胀系数等热物性，并且与实验数据呈现良好的一致性. Bandrés等人[6]研究了离子液体在不同的压力、温度下的黏度，采用分子动力学模拟方法，从分子角度分析了静电、分子大小、阴阳离子间相互作用强度对液体黏度的影响. Doherty等人[7]模拟了5种DES的密度、黏度和表面张力，通过优化模拟体系的OPLS力场参数，取得了与实验吻合的结果，但优化后的力场仅限于这5种特定DES体系的模拟，无法应用于其他低共熔溶剂的模拟. 从现有文献来看，利用模拟技术研究低共熔溶剂的物性和结构的报道仍比较缺乏，并且其微观相互作用机制与物性之间的关系尚未完全弄清楚.

本文采用分子动力学模拟的方法，研究氯化胆碱类低共熔溶剂在不同温度下的密度和黏度，通过计算径向分布函数，从微观上研究其行为变化对低共熔溶剂物性产生的影响，指导合成特定性质的低共熔溶剂，为其科学研究和工业应用提供支持.

2 模型与方法

2.1 模拟体系

本文选用了3种不同的DES体系，每种体系名称、组成成分和摩尔比如表1所示. 其中各种氢键受体(HBA)、氢键供体(HBD)的结构如图1所示.

表1 低共熔溶剂体系组成

(a) 氯化胆碱(ChCl)

(b) 乙二醇(ETG)

(c) 1,3-丙二醇(1,3-PRO)

(d)1,4-丁二醇(1,4-BDL)

当氯化胆碱(ChCl)和不同的氢键供体(HBD)配对进行分子动力学模拟时，为了方便后续分析不同原子间的相互作用机理，将氯化胆碱中的C原子用Cc标识，H原子用Hc标识，O原子用Oc标识，用来区分HBD中的C、H、O原子.

2.2 模拟方法

本文使用的软件为GROMACS 2016.4[8]版本，选用的力场为GAFF (Generation Amber Force Field)力场[9]，这种力场处理各种有机小分子具有优势[10]. 由于GROMACS程序没有自带的GAFF力场可用，因此需要使用AmberTools[11]产生模拟所需要的力场参数，再通过ACPYPE[12]将其转换为GROMACS形式的可用参数，方可进行模拟. 采用Lennard-Jones 12-6(LJ)势函数计算力场中的范德华相互作用，其表达式为：

(1)

LJ势中的参数采用Lorentz-Berthelot几何混合规则来计算：

(2)

(3)

其中，rij代表原子i和j之间的距离；εij，σij为不同原子间的势阱深度和尺寸参数[13].

由于低共熔溶剂是类离子液体，体系中存在游离的氯离子以及部分带电荷的基团，不同的方法计算获得的电荷对模拟结果有很大影响，因此电荷的计算方法显得尤为重要. 本文采用了限制性静电势拟合(restrained electrostatic potential, RESP)[14]的方法获取体系的电荷. 这种方法解决了传统的原子电荷计算方法存在的缺陷，通过对每个原子加以限制作用避免因原子距离过远、对构象依赖大等问题造成的拟合不准确，计算的电荷质量相对较高[15]. 文中采用Multiwfn 3.7[16]程序计算RESP电荷.分子动力学模拟类离子液体时，其静电势通常被高估[17]，可用控制静电相互作用的方法将计算获得的电荷减小5%～25%[18]. 对于DES体系，选择降低10%，即采用0.9的比例因子来修正计算的RESP电荷获得了较好的效果.

首先构建模拟体系所需分子的结构图(如图1所示)，然后使用Packmol[19]创建并填充立方体盒子，尺寸为6 nm×6 nm×6 nm，并按照给定比例将不同的分子随机插入到立方盒子中，填充的分子数如表1所示. 在填充氯化胆碱(ChCl)时，需要分别构建氯化胆碱阳离子([Cho]+)和氯离子(Cl-)，将其各自独立地插入到模拟盒子中，因为在低共熔溶剂中，氯离子是以游离状态存在的，将其独立地插入盒子中，有利于低共熔溶剂各组分的充分混合，使模拟状态更接近实际状态. 对模拟体系的x、y、z方向均设置周期性边界条件.首先进行能量最小化，防止构建的体系原子分布不合理.采用的方法为最速下降法，初始步长为0.01 nm，设定最大力小于1.0 kJ·mol-1·nm-1时认为最小化达到收敛. 选择cut-off方式计算范德华相互作用，PME(Particle Mesh Ewald)计算原子静电相互作用，库伦截断距离和范德华截断距离均为1.0 nm. 然后进行系统平衡，压力为0.1 MPa，平衡期间压力变化比较大，因此采用Berendsen方法控压，V-rescale方法进行控温，时间常数均为1.0 ps. 积分步长为0.002 ps，总平衡时间为5 ns，在平衡阶段的前100 ps对系统做了退火处理，按照指定速度缓慢升温，防止温度变化时原子速度变化过快，体系不稳定. 最后进行模拟计算，时间为10 ns，每500步输出一次轨迹文件，控压方法改为parrinello-rahman方式，适用于平衡后的体系控压. 用lincs算法对含有H原子的键施加约束. 密度和黏度的模拟都是在NPT系综下进行的，黏度的计算采用的是周期性扰动法[20].

3 结果与讨论

3.1 密 度

表2给出了293.15～353.15 K温度范围内模拟获得的低共熔溶剂的密度.每个温度点重复模拟三次后取平均值，并和文献值进行比较. 图2所示为DES密度和温度的关系. 当氯化胆碱和相同的HBD组合时，HBD越多，形成的DES体系密度越小，但不同摩尔比下的密度差值均保持在6 kg·m-3以内，随摩尔比变化程度不大. 密度与温度基本呈线性关系，并且每个数据点均落在各自的拟合线上，因此DES的密度与温度的关系可以用公式(4)进行描述[21]：

ρ=AT+B

(4)

其中，T为绝对温度(K)；A、B均为拟合参数. 每种DES的拟合参数如表3所示.

3.2 黏 度

通过模拟出的每种DES体系的黏度与文献实验值的误差如表2所示，黏度随温度的变化如图3所示，均随温度升高明显降低. 黏度大小为Etgline<13proline (1∶4)<13proline (1∶3)<14bdline (1∶4)<14bdline (1∶3). 其中14bdline (1∶3)体系的黏度随温度变化程度最大，温度敏感性最高. 且HBD的占比越多，形成的DES体系黏度越小. DES黏度随温度的显著变化可以用阿伦尼乌斯方程 (Arrhenius equation)来描述[22]：

(5)

其中，T为绝对温度(K);η0为常数；Eη为表观活化能(kJ·mol-1)；R为摩尔气体常数(kJ·mol-1·K-1).不同DES体系的参数如表3所示，相关系数均在0.99以上，说明此方程可以很好地描述DES黏度与温度的关系.

表2 氯化胆碱+乙二醇/1，3-丙二醇/1，4-丁二醇体系的密度、黏度的模拟值与实验值及其相对偏差

表3 不同DES体系关于公式(4)和公式(5)的参数

图2 不同DES的密度与温度的关系Fig.2 Relationship between density and temperature of different DES

图3 不同DES的黏度与温度的关系Fig.3 Relationship between viscosity and temperature of different DES

3.3 结构分析

径向分布函数是描述分子结构常用的方法. 本文采用GAFF力场模拟出的密度和黏度与实验值吻合较好，证明了所使用的力场参数以及电荷计算方法适用于所选的低共熔溶剂，可以采用径向分布函数研究低共熔溶剂体系的内部结构.

径向分布函数是指液体中距离某一参考原子A周围r处原子B的密度，它反映了系统中短程有序，长程无序的特点. 当参考原子为中心原子A时，距离原子A中心r到r+dr范围内的分子数为dN，则径向分布函数g(r)[13]表示为：

ρg(r)4πr2=dN

(6)

(7)

其中ρ为系统中液体密度；N为系统中分子数.因此，径向分布函数也可以理解为系统的局部密度和平均密度的比值，表示参考原子A周围原子B出现的概率[11].

3.3.1 氢键主要作用位点 为了考察DES体系中氢键的主要存在位点，我们对其径向分布函数进行了分析. 本文所研究的体系均为氯化胆碱+二元醇，且二元醇的两个羟基均位于两端. 计算发现不同体系中同类基团之间的径向分布函数趋势相似. 为了简洁起见，我们以Etgline体系为例进行主要分析. 由图4a所示，当参考原子为Cl时，比较明显的峰出现在Hc21和H8原子上，分别为[Cho]+的羟基氢原子Hc21和HBD的羟基氢原子H8，距离分别为2.26 和2.30 Å，均小于Cl和H之间的原子范德华半径之和2.86 Å[26]. GROMACS中氢键判断的标准为：原子间距小于3.5 Å，且峰形比较尖锐. 因此，在这两处氢原子存在较强的、相对稳定的相互作用，可判定为氢键作用. 相比之下，其余两种甲基上的氢原子：Hc2([Cho]+)、H2(HBD)与阴离子间的最短距离分别为2.9和3.56 Å，大于原子范德华半径之和，且峰值明显较低，不能判定向阴离子提供氢键，以范德华作用形式存在. 综上分析，由Etgline体系的径向分布函数图可知，以Cl为参考原子时，-OH上的氢原子多在2.3 Å左右处聚集，原子之间跨越了范德华表面而形成了强氢键相互作用；其余氢原子多在3～6 Å处以范德华相互作用聚集，相互作用相对较弱；当距离大于7 Å时，体系中的各原子呈无序状态排列.

由图4b所示，以乙二醇的羟基氧原子O7为参考原子时，考察[Cho]+上的羟基氢原子(Hc21)和乙二醇的羟基氢原子(H8、H10)的相互作用. O7-H8之间出现了尖锐的RDF峰，主要是由于这两个原子直接以键连方式存在，相互作用明显最强，由于坐标范围有限，而O7-H8之间由于键连关系形成的峰高且窄，因此在图中的显示如同直线. O7与H10、Hc21的峰值相对较低，虽然O的电负性(3.44)比Cl(3.16)大，但在HBD中，O和H相连以羟基方式存在，受到H原子影响，-OH整体的电负性大大降低且小于Cl，因此其聚集吸引H原子的能力也不如Cl，相互作用明显减弱. 但O7与H10、Hc21的最短接触距离分别为1.84和1.80 Å，小于氢键作用范围，因此可判定为氢键供体的作用位点. O7和Cl为参考原子时分别呈现多个峰和单峰，原因是O7同时存在分子内和分子间的作用力，而Cl是游离状态，只存在分子间作用力.

图5是Etgline中各部分分子质心间(COM)的径向分布函数. 从图5中可以看出，Cl与[Cho]+之间相互作用最强，其次是Cl和HBD. 根据图5a、5b可知，在6 Å以内，不同分子(Cl-[Cho]+、Cl-HBD、HBD-HBD)之间的相互作用占主导地位，所有分子间(包括Cl-Cl、[Cho]+-[Cho]+、HBD-HBD)的最近接触距离均保持在8 Å以内.

综上，氯化胆碱和乙二醇组成的低共熔溶剂Etgline体系中，氢键网络主要是由Cl原子和[Cho]+的羟基氢(Hc21)、HBD的羟基氢(H8)原子形成的，且Cl与Hc21 ([Cho]+)的相互作用最强.

3.3.2 温度对DES体系行为的影响 利用径向分布函数，本文还研究了温度对氢键体系的影响. 由于氢键主要存在于阴离子和[Cho]+、HBD上的羟基氢原子，因此选择计算的径向分布函数为以下两组：Cl-H([Cho]+)、Cl-H (HBD).

图6、图7和图8所示为由氢键受体(HBA)为氯化胆碱，氢键供体(HBD)分别为乙二醇、1，3-丙二醇、1，4-丁二醇所形成的低共熔溶剂，在温度为273.15、313.15和353.15 K，压力为0.1 MPa时的Cl-Hc21 ([Cho]+)、Cl-H (HBD)原子间的径向分布函数. 从图中可以看到，当温度升高时，Etgline、1,3proline、1,4bdline三种体系中，阴离子与不同氢原子之间呈现的径向分布函数峰值均有所下降，表明温度影响体系的氢键网络，导致分子间的相互作用强度减弱; 但不同温度下峰值对应的原子间距基本没有变化，表明氢键网络的结构基本未发生变化.

表4 同DES中Cl-Hc21, Cl-H在不同温度下的峰高和最近距离

从表4可以得出以下结论：(1) 每个体系中不同温度下对应的Cl-Hc21 ([Cho]+)形成的氢键键长基本稳定在2.26 Å，Cl-H (HBD)键长基本稳定在2.3 Å；(2) 随着温度升高，Cl-Hc21和Cl-H (HBD)的峰值均逐渐减小，即Cl周围存在的氯化胆碱阳离子上的羟基氢原子减少，可推断升高温度，可能会使Cl和H之间的相互作用减弱，对氢键的形成有破坏；(3) 温度升高破坏分子间相互作用力，体系中部分分子脱离了强氢键作用的束缚，更容易扩散，造成体系黏度降低.

3.3.3 烷基链长度对DES体系行为的影响 为了研究作为氢键供体的烷基链长度对体系结构的影响，计算的Cl-H (HBD)、Cl-H ([Cho]+)原子间的径向分布函数如图9所示. 为了排除HBD相对数量的影响，计算时采用的HBA∶HBD摩尔比均为1∶4.

从图9可以看出，随着氢键供体烷基链长度的增加，Cl原子与氯化胆碱阳离子、氢键供体上的羟基氢原子(Cl-H([Cho]+)、Cl-H(HBD))之间的相互作用均增强. 根据之前模拟的数据可知，黏度大小为ChCl+乙二醇

4 结 论

本文采用了分子动力学模拟的方法，研究了氯化胆碱+二元醇类低共熔溶剂的性质和结构特性，主要结论有：(1) 对于氯化胆碱与乙二醇/1，3-丙二醇/1，4-丁二醇体系，采用GAFF力场，0.9RESP电荷模拟的密度、黏度值与文献中实验值的平均绝对偏差分别为0.30%和6.19%，可以很好地重现实验数据. 混合形成的低共熔溶剂的密度和黏度均随温度升高而降低，密度随温度呈线性变化，黏度和温度的关系符合阿伦尼乌斯定律，密度和黏度拟合方程的计算值与实验值的平均绝对偏差分别为0.27%和7.99%，具有良好的预测性. (2) 低共熔溶剂的低熔点与组分之间形成广泛的氢键有关，而氢键体系主要是由Cl原子与-OH(包括氢键供体和氢键受体)上的氢形成的. (3) 根据RDF可知：随着温度升高，热运动剧烈，分子间相互作用力减弱，对氢键有一定的破坏，导致黏度的降低；随着HBD烷基链长度的增加，同类原子间相互作用力增强，黏度增大.