聶文潔， 王劍飛， 趙貫甲， 尹建國， 馬素霞

(1. 太原理工大學(xué)電氣與動力工程學(xué)院， 太原 030024；2. 太原理工大學(xué)循環(huán)流化床高效清潔與利用山西省重點實驗室， 太原 030024；3.太原鍋爐集團有限公司， 太原 030024)

1 引 言

低共熔溶劑(deep eutectic solvents)是一種新型的綠色溶劑，不僅具備離子液體的蒸汽壓低、溶解性好、導(dǎo)電性強、化學(xué)穩(wěn)定性好等諸多優(yōu)勢，還克服了離子液體制備復(fù)雜、提純困難、微毒、價格昂貴、不易降解等缺點[1]. 2003年，Abbott等[2]首次將氯化膽堿(choline chloride，ChCl)和尿素(urea)在1∶2的摩爾比下加熱攪拌混合，冷卻后形成了一種熔點很低的混合液體，即低共熔溶劑.由于低共熔溶劑制備簡單、原料綠色廉價，且兼具離子液體良好的化學(xué)性質(zhì)，引起了廣大研究者的重視. 其在酸性氣體吸收、分離過程、化學(xué)反應(yīng)、生物催化和電化學(xué)等領(lǐng)域顯示出良好的應(yīng)用前景[3, 4].

低共熔溶劑至少由兩種成分組成：一種稱為氫鍵受體(HBA)，一種稱為氫鍵供體(HBD).氫鍵受體通常有季銨鹽(如氯化膽堿)和兩性離子(如甜菜堿)等；氫鍵供體主要有尿素、醇類(如乙二醇、甘油、丙二醇、丁二醇)、酸類(乙酰丙酸、丙二酸、草酸)和苯酚等[1]. 由于每種HBA和HBD都能隨意組合，因此理論上可形成無限多種低共熔溶劑. 其組成不同，結(jié)構(gòu)和性質(zhì)也不同. 氫鍵的形成是導(dǎo)致低共熔溶劑熔點降低和宏觀熱力學(xué)性質(zhì)變化的重要因素. 分子動力學(xué)模擬提供了一種可以計算流體宏觀特性和微觀結(jié)構(gòu)的方法，是研究這類體系經(jīng)濟、可靠的手段. Logotheti等人[5]采用分子動力學(xué)模擬計算了298.15～333.15 K溫度下甲基咪唑類離子液體的密度、擴散系數(shù)和熱膨脹系數(shù)等熱物性，并且與實驗數(shù)據(jù)呈現(xiàn)良好的一致性. Bandrés等人[6]研究了離子液體在不同的壓力、溫度下的黏度，采用分子動力學(xué)模擬方法，從分子角度分析了靜電、分子大小、陰陽離子間相互作用強度對液體黏度的影響. Doherty等人[7]模擬了5種DES的密度、黏度和表面張力，通過優(yōu)化模擬體系的OPLS力場參數(shù)，取得了與實驗吻合的結(jié)果，但優(yōu)化后的力場僅限于這5種特定DES體系的模擬，無法應(yīng)用于其他低共熔溶劑的模擬. 從現(xiàn)有文獻來看，利用模擬技術(shù)研究低共熔溶劑的物性和結(jié)構(gòu)的報道仍比較缺乏，并且其微觀相互作用機制與物性之間的關(guān)系尚未完全弄清楚.

本文采用分子動力學(xué)模擬的方法，研究氯化膽堿類低共熔溶劑在不同溫度下的密度和黏度，通過計算徑向分布函數(shù)，從微觀上研究其行為變化對低共熔溶劑物性產(chǎn)生的影響，指導(dǎo)合成特定性質(zhì)的低共熔溶劑，為其科學(xué)研究和工業(yè)應(yīng)用提供支持.

2 模型與方法

2.1 模擬體系

本文選用了3種不同的DES體系，每種體系名稱、組成成分和摩爾比如表1所示. 其中各種氫鍵受體(HBA)、氫鍵供體(HBD)的結(jié)構(gòu)如圖1所示.

表1 低共熔溶劑體系組成

(a) 氯化膽堿(ChCl)

(b) 乙二醇(ETG)

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

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

當(dāng)氯化膽堿(ChCl)和不同的氫鍵供體(HBD)配對進行分子動力學(xué)模擬時，為了方便后續(xù)分析不同原子間的相互作用機理，將氯化膽堿中的C原子用Cc標識，H原子用Hc標識，O原子用Oc標識，用來區(qū)分HBD中的C、H、O原子.

2.2 模擬方法

本文使用的軟件為GROMACS 2016.4[8]版本，選用的力場為GAFF (Generation Amber Force Field)力場[9]，這種力場處理各種有機小分子具有優(yōu)勢[10]. 由于GROMACS程序沒有自帶的GAFF力場可用，因此需要使用AmberTools[11]產(chǎn)生模擬所需要的力場參數(shù)，再通過ACPYPE[12]將其轉(zhuǎn)換為GROMACS形式的可用參數(shù)，方可進行模擬. 采用Lennard-Jones 12-6(LJ)勢函數(shù)計算力場中的范德華相互作用，其表達式為：

(1)

LJ勢中的參數(shù)采用Lorentz-Berthelot幾何混合規(guī)則來計算：

(2)

(3)

其中，rij代表原子i和j之間的距離；εij，σij為不同原子間的勢阱深度和尺寸參數(shù)[13].

由于低共熔溶劑是類離子液體，體系中存在游離的氯離子以及部分帶電荷的基團，不同的方法計算獲得的電荷對模擬結(jié)果有很大影響，因此電荷的計算方法顯得尤為重要. 本文采用了限制性靜電勢擬合(restrained electrostatic potential, RESP)[14]的方法獲取體系的電荷. 這種方法解決了傳統(tǒng)的原子電荷計算方法存在的缺陷，通過對每個原子加以限制作用避免因原子距離過遠、對構(gòu)象依賴大等問題造成的擬合不準確，計算的電荷質(zhì)量相對較高[15]. 文中采用Multiwfn 3.7[16]程序計算RESP電荷.分子動力學(xué)模擬類離子液體時，其靜電勢通常被高估[17]，可用控制靜電相互作用的方法將計算獲得的電荷減小5%～25%[18]. 對于DES體系，選擇降低10%，即采用0.9的比例因子來修正計算的RESP電荷獲得了較好的效果.

首先構(gòu)建模擬體系所需分子的結(jié)構(gòu)圖(如圖1所示)，然后使用Packmol[19]創(chuàng)建并填充立方體盒子，尺寸為6 nm×6 nm×6 nm，并按照給定比例將不同的分子隨機插入到立方盒子中，填充的分子數(shù)如表1所示. 在填充氯化膽堿(ChCl)時，需要分別構(gòu)建氯化膽堿陽離子([Cho]+)和氯離子(Cl-)，將其各自獨立地插入到模擬盒子中，因為在低共熔溶劑中，氯離子是以游離狀態(tài)存在的，將其獨立地插入盒子中，有利于低共熔溶劑各組分的充分混合，使模擬狀態(tài)更接近實際狀態(tài). 對模擬體系的x、y、z方向均設(shè)置周期性邊界條件.首先進行能量最小化，防止構(gòu)建的體系原子分布不合理.采用的方法為最速下降法，初始步長為0.01 nm，設(shè)定最大力小于1.0 kJ·mol-1·nm-1時認為最小化達到收斂. 選擇cut-off方式計算范德華相互作用，PME(Particle Mesh Ewald)計算原子靜電相互作用，庫倫截斷距離和范德華截斷距離均為1.0 nm. 然后進行系統(tǒng)平衡，壓力為0.1 MPa，平衡期間壓力變化比較大，因此采用Berendsen方法控壓，V-rescale方法進行控溫，時間常數(shù)均為1.0 ps. 積分步長為0.002 ps，總平衡時間為5 ns，在平衡階段的前100 ps對系統(tǒng)做了退火處理，按照指定速度緩慢升溫，防止溫度變化時原子速度變化過快，體系不穩(wěn)定. 最后進行模擬計算，時間為10 ns，每500步輸出一次軌跡文件，控壓方法改為parrinello-rahman方式，適用于平衡后的體系控壓. 用lincs算法對含有H原子的鍵施加約束. 密度和黏度的模擬都是在NPT系綜下進行的，黏度的計算采用的是周期性擾動法[20].

3 結(jié)果與討論

3.1 密 度

表2給出了293.15～353.15 K溫度范圍內(nèi)模擬獲得的低共熔溶劑的密度.每個溫度點重復(fù)模擬三次后取平均值，并和文獻值進行比較. 圖2所示為DES密度和溫度的關(guān)系. 當(dāng)氯化膽堿和相同的HBD組合時，HBD越多，形成的DES體系密度越小，但不同摩爾比下的密度差值均保持在6 kg·m-3以內(nèi)，隨摩爾比變化程度不大. 密度與溫度基本呈線性關(guān)系，并且每個數(shù)據(jù)點均落在各自的擬合線上，因此DES的密度與溫度的關(guān)系可以用公式(4)進行描述[21]：

ρ=AT+B

(4)

其中，T為絕對溫度(K)；A、B均為擬合參數(shù). 每種DES的擬合參數(shù)如表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為常數(shù)；Eη為表觀活化能(kJ·mol-1)；R為摩爾氣體常數(shù)(kJ·mol-1·K-1).不同DES體系的參數(shù)如表3所示，相關(guān)系數(shù)均在0.99以上，說明此方程可以很好地描述DES黏度與溫度的關(guān)系.

表2 氯化膽堿+乙二醇/1，3-丙二醇/1，4-丁二醇體系的密度、黏度的模擬值與實驗值及其相對偏差

表3 不同DES體系關(guān)于公式(4)和公式(5)的參數(shù)

圖2 不同DES的密度與溫度的關(guān)系Fig.2 Relationship between density and temperature of different DES

圖3 不同DES的黏度與溫度的關(guān)系Fig.3 Relationship between viscosity and temperature of different DES

3.3 結(jié)構(gòu)分析

徑向分布函數(shù)是描述分子結(jié)構(gòu)常用的方法. 本文采用GAFF力場模擬出的密度和黏度與實驗值吻合較好，證明了所使用的力場參數(shù)以及電荷計算方法適用于所選的低共熔溶劑，可以采用徑向分布函數(shù)研究低共熔溶劑體系的內(nèi)部結(jié)構(gòu).

徑向分布函數(shù)是指液體中距離某一參考原子A周圍r處原子B的密度，它反映了系統(tǒng)中短程有序，長程無序的特點. 當(dāng)參考原子為中心原子A時，距離原子A中心r到r+dr范圍內(nèi)的分子數(shù)為dN，則徑向分布函數(shù)g(r)[13]表示為：

ρg(r)4πr2=dN

(6)

(7)

其中ρ為系統(tǒng)中液體密度；N為系統(tǒng)中分子數(shù).因此，徑向分布函數(shù)也可以理解為系統(tǒng)的局部密度和平均密度的比值，表示參考原子A周圍原子B出現(xiàn)的概率[11].

3.3.1 氫鍵主要作用位點 為了考察DES體系中氫鍵的主要存在位點，我們對其徑向分布函數(shù)進行了分析. 本文所研究的體系均為氯化膽堿+二元醇，且二元醇的兩個羥基均位于兩端. 計算發(fā)現(xiàn)不同體系中同類基團之間的徑向分布函數(shù)趨勢相似. 為了簡潔起見，我們以Etgline體系為例進行主要分析. 由圖4a所示，當(dāng)參考原子為Cl時，比較明顯的峰出現(xiàn)在Hc21和H8原子上，分別為[Cho]+的羥基氫原子Hc21和HBD的羥基氫原子H8，距離分別為2.26 和2.30 ?，均小于Cl和H之間的原子范德華半徑之和2.86 ?[26]. GROMACS中氫鍵判斷的標準為：原子間距小于3.5 ?，且峰形比較尖銳. 因此，在這兩處氫原子存在較強的、相對穩(wěn)定的相互作用，可判定為氫鍵作用. 相比之下，其余兩種甲基上的氫原子：Hc2([Cho]+)、H2(HBD)與陰離子間的最短距離分別為2.9和3.56 ?，大于原子范德華半徑之和，且峰值明顯較低，不能判定向陰離子提供氫鍵，以范德華作用形式存在. 綜上分析，由Etgline體系的徑向分布函數(shù)圖可知，以Cl為參考原子時，-OH上的氫原子多在2.3 ?左右處聚集，原子之間跨越了范德華表面而形成了強氫鍵相互作用；其余氫原子多在3～6 ?處以范德華相互作用聚集，相互作用相對較弱；當(dāng)距離大于7 ?時，體系中的各原子呈無序狀態(tài)排列.

由圖4b所示，以乙二醇的羥基氧原子O7為參考原子時，考察[Cho]+上的羥基氫原子(Hc21)和乙二醇的羥基氫原子(H8、H10)的相互作用. O7-H8之間出現(xiàn)了尖銳的RDF峰，主要是由于這兩個原子直接以鍵連方式存在，相互作用明顯最強，由于坐標范圍有限，而O7-H8之間由于鍵連關(guān)系形成的峰高且窄，因此在圖中的顯示如同直線. 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為參考原子時分別呈現(xiàn)多個峰和單峰，原因是O7同時存在分子內(nèi)和分子間的作用力，而Cl是游離狀態(tài)，只存在分子間作用力.

圖5是Etgline中各部分分子質(zhì)心間(COM)的徑向分布函數(shù). 從圖5中可以看出，Cl與[Cho]+之間相互作用最強，其次是Cl和HBD. 根據(jù)圖5a、5b可知，在6 ?以內(nèi)，不同分子(Cl-[Cho]+、Cl-HBD、HBD-HBD)之間的相互作用占主導(dǎo)地位，所有分子間(包括Cl-Cl、[Cho]+-[Cho]+、HBD-HBD)的最近接觸距離均保持在8 ?以內(nèi).

綜上，氯化膽堿和乙二醇組成的低共熔溶劑Etgline體系中，氫鍵網(wǎng)絡(luò)主要是由Cl原子和[Cho]+的羥基氫(Hc21)、HBD的羥基氫(H8)原子形成的，且Cl與Hc21 ([Cho]+)的相互作用最強.

3.3.2 溫度對DES體系行為的影響 利用徑向分布函數(shù)，本文還研究了溫度對氫鍵體系的影響. 由于氫鍵主要存在于陰離子和[Cho]+、HBD上的羥基氫原子，因此選擇計算的徑向分布函數(shù)為以下兩組：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)原子間的徑向分布函數(shù). 從圖中可以看到，當(dāng)溫度升高時，Etgline、1,3proline、1,4bdline三種體系中，陰離子與不同氫原子之間呈現(xiàn)的徑向分布函數(shù)峰值均有所下降，表明溫度影響體系的氫鍵網(wǎng)絡(luò)，導(dǎo)致分子間的相互作用強度減弱; 但不同溫度下峰值對應(yīng)的原子間距基本沒有變化，表明氫鍵網(wǎng)絡(luò)的結(jié)構(gòu)基本未發(fā)生變化.

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

從表4可以得出以下結(jié)論：(1) 每個體系中不同溫度下對應(yīng)的Cl-Hc21 ([Cho]+)形成的氫鍵鍵長基本穩(wěn)定在2.26 ?，Cl-H (HBD)鍵長基本穩(wěn)定在2.3 ?；(2) 隨著溫度升高，Cl-Hc21和Cl-H (HBD)的峰值均逐漸減小，即Cl周圍存在的氯化膽堿陽離子上的羥基氫原子減少，可推斷升高溫度，可能會使Cl和H之間的相互作用減弱，對氫鍵的形成有破壞；(3) 溫度升高破壞分子間相互作用力，體系中部分分子脫離了強氫鍵作用的束縛，更容易擴散，造成體系黏度降低.

3.3.3 烷基鏈長度對DES體系行為的影響 為了研究作為氫鍵供體的烷基鏈長度對體系結(jié)構(gòu)的影響，計算的Cl-H (HBD)、Cl-H ([Cho]+)原子間的徑向分布函數(shù)如圖9所示. 為了排除HBD相對數(shù)量的影響，計算時采用的HBA∶HBD摩爾比均為1∶4.

從圖9可以看出，隨著氫鍵供體烷基鏈長度的增加，Cl原子與氯化膽堿陽離子、氫鍵供體上的羥基氫原子(Cl-H([Cho]+)、Cl-H(HBD))之間的相互作用均增強. 根據(jù)之前模擬的數(shù)據(jù)可知，黏度大小為ChCl+乙二醇

4 結(jié) 論

本文采用了分子動力學(xué)模擬的方法，研究了氯化膽堿+二元醇類低共熔溶劑的性質(zhì)和結(jié)構(gòu)特性，主要結(jié)論有：(1) 對于氯化膽堿與乙二醇/1，3-丙二醇/1，4-丁二醇體系，采用GAFF力場，0.9RESP電荷模擬的密度、黏度值與文獻中實驗值的平均絕對偏差分別為0.30%和6.19%，可以很好地重現(xiàn)實驗數(shù)據(jù). 混合形成的低共熔溶劑的密度和黏度均隨溫度升高而降低，密度隨溫度呈線性變化，黏度和溫度的關(guān)系符合阿倫尼烏斯定律，密度和黏度擬合方程的計算值與實驗值的平均絕對偏差分別為0.27%和7.99%，具有良好的預(yù)測性. (2) 低共熔溶劑的低熔點與組分之間形成廣泛的氫鍵有關(guān)，而氫鍵體系主要是由Cl原子與-OH(包括氫鍵供體和氫鍵受體)上的氫形成的. (3) 根據(jù)RDF可知：隨著溫度升高，熱運動劇烈，分子間相互作用力減弱，對氫鍵有一定的破壞，導(dǎo)致黏度的降低；隨著HBD烷基鏈長度的增加，同類原子間相互作用力增強，黏度增大.