氯化膽堿類(lèi)低共熔溶劑結(jié)構(gòu)與性質(zhì)的分子動(dòng)力學(xué)研究

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

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

1 引 言

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

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

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

2 模型與方法

2.1 模擬體系

本文選用了3種不同的DES體系，每種體系名稱(chēng)、組成成分和摩爾比如表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)配對(duì)進(jìn)行分子動(dòng)力學(xué)模擬時(shí)，為了方便后續(xù)分析不同原子間的相互作用機(jī)理，將氯化膽堿中的C原子用Cc標(biāo)識(shí)，H原子用Hc標(biāo)識(shí)，O原子用Oc標(biāo)識(shí)，用來(lái)區(qū)分HBD中的C、H、O原子.

2.2 模擬方法

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

(1)

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

(2)

(3)

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

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

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

3 結(jié)果與討論

3.1 密 度

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

ρ=AT+B

(4)

其中，T為絕對(duì)溫度(K)；A、B均為擬合參數(shù). 每種DES的擬合參數(shù)如表3所示.

3.2 黏 度

通過(guò)模擬出的每種DES體系的黏度與文獻(xiàn)實(shí)驗(yàn)值的誤差如表2所示，黏度隨溫度的變化如圖3所示，均隨溫度升高明顯降低. 黏度大小為Etgline<13proline (1∶4)<13proline (1∶3)<14bdline (1∶4)<14bdline (1∶3). 其中14bdline (1∶3)體系的黏度隨溫度變化程度最大，溫度敏感性最高. 且HBD的占比越多，形成的DES體系黏度越小. DES黏度隨溫度的顯著變化可以用阿倫尼烏斯方程 (Arrhenius equation)來(lái)描述[22]：

(5)

其中，T為絕對(duì)溫度(K);η0為常數(shù)；Eη為表觀(guān)活化能(kJ·mol-1)；R為摩爾氣體常數(shù)(kJ·mol-1·K-1).不同DES體系的參數(shù)如表3所示，相關(guān)系數(shù)均在0.99以上，說(shuō)明此方程可以很好地描述DES黏度與溫度的關(guān)系.

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

表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力場(chǎng)模擬出的密度和黏度與實(shí)驗(yàn)值吻合較好，證明了所使用的力場(chǎng)參數(shù)以及電荷計(jì)算方法適用于所選的低共熔溶劑，可以采用徑向分布函數(shù)研究低共熔溶劑體系的內(nèi)部結(jié)構(gòu).

徑向分布函數(shù)是指液體中距離某一參考原子A周?chē)鷕處原子B的密度，它反映了系統(tǒng)中短程有序，長(zhǎng)程無(wú)序的特點(diǎn). 當(dāng)參考原子為中心原子A時(shí)，距離原子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周?chē)覤出現(xiàn)的概率[11].

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

由圖4b所示，以乙二醇的羥基氧原子O7為參考原子時(shí)，考察[Cho]+上的羥基氫原子(Hc21)和乙二醇的羥基氫原子(H8、H10)的相互作用. O7-H8之間出現(xiàn)了尖銳的RDF峰，主要是由于這兩個(gè)原子直接以鍵連方式存在，相互作用明顯最強(qiáng)，由于坐標(biāo)范圍有限，而O7-H8之間由于鍵連關(guān)系形成的峰高且窄，因此在圖中的顯示如同直線(xiàn). O7與H10、Hc21的峰值相對(duì)較低，雖然O的電負(fù)性(3.44)比Cl(3.16)大，但在HBD中，O和H相連以羥基方式存在，受到H原子影響，-OH整體的電負(fù)性大大降低且小于Cl，因此其聚集吸引H原子的能力也不如Cl，相互作用明顯減弱. 但O7與H10、Hc21的最短接觸距離分別為1.84和1.80 ?，小于氫鍵作用范圍，因此可判定為氫鍵供體的作用位點(diǎn). O7和Cl為參考原子時(shí)分別呈現(xiàn)多個(gè)峰和單峰，原因是O7同時(shí)存在分子內(nèi)和分子間的作用力，而Cl是游離狀態(tài)，只存在分子間作用力.

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

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

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

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

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

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

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

4 結(jié) 論

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