水/石墨烯界面離子吸附的分子動(dòng)力學(xué)模擬: 力場參數(shù)優(yōu)化與吸附機(jī)制

廖首維， 劉炎昌， 石澤南， 趙道輝， 魏嫣瑩，3， 李理波，3

（1. 華南理工大學(xué)化學(xué)與化工學(xué)院, 廣東省綠色化學(xué)產(chǎn)品技術(shù)重點(diǎn)實(shí)驗(yàn)室, 廣州 510640；2. 湖北大學(xué)化學(xué)與化工學(xué)院, 武漢 430062； 3. 華南理工大學(xué)制漿造紙工程國家重點(diǎn)實(shí)驗(yàn)室, 廣州 510640）

石墨烯、 二硫化鉬和Mxene等二維材料［1～3］構(gòu)建的納米孔道在離子輸運(yùn)與篩分、 水輸運(yùn)等方面表現(xiàn)出大量的新穎現(xiàn)象和優(yōu)異性能［4～8］， 因而在許多領(lǐng)域（如海水淡化［5，9～13］、 滲透發(fā)電［14～17］、 超級(jí)電容器［18～20］、 納米元器件［21，22］等）具有廣闊的應(yīng)用前景. 與之相關(guān)的研究往往涉及鹽溶液-二維材料體系， 時(shí)空尺度很微小， 通常的實(shí)驗(yàn)技術(shù)難以靈敏檢測(cè)［23］.

量子化學(xué)需要的計(jì)算資源太多， 通常只能模擬數(shù)百原子、 數(shù)百皮秒的時(shí)空尺度［24］； 傳統(tǒng)連續(xù)介質(zhì)模型精度不夠高， 很難準(zhǔn)確描述1 nm以下的尺度［25，26］： 它們都很難用于鹽溶液-二維材料體系的研究.經(jīng)典分子動(dòng)力學(xué)（Molecular dynamics， MD）模擬因其處理的時(shí)空尺度合適（顯著大于量子化學(xué)計(jì)算）、 模擬精度較高（顯著高于連續(xù)介質(zhì)模型）， 已被廣泛應(yīng)用于各種鹽溶液-二維材料體系［20，27～31］的研究， 為實(shí)驗(yàn)研究提供越來越重要的指導(dǎo)和補(bǔ)充作用.

目前， 用MD模擬研究二維材料納米孔道里鹽溶液面臨的主要問題有： 首先， 石墨烯等金屬、 半導(dǎo)體二維材料與離子之間存在強(qiáng)烈的離子-π相互作用［32～37］. 最近的模擬［38］和實(shí)驗(yàn)研究［9，37］都表明， 離子-π相互作用對(duì)離子在石墨烯納米孔道里的吸附與輸運(yùn)很重要， 然而當(dāng)前的固定電荷力場（Fixed charge force field）設(shè)定石墨烯C原子電荷為0（因?yàn)楦鰿原子化學(xué)環(huán)境相同）， 無法體現(xiàn)離子-π相互作用引起的電荷變化， 導(dǎo)致無法模擬離子在石墨烯表面的吸附. 從頭算分子動(dòng)力學(xué)（AIMD）［39］可模擬的時(shí)空尺度太?。?可極化力場MD模擬［35］的運(yùn)算速度也很慢， 且該方法還在繼續(xù)發(fā)展中. 而MD模擬力場中不同原子之間的Lennard-Jones（LJ）作用參數(shù)往往使用混合規(guī)則（如Lorentz-Berthelot 或LB 規(guī)則）計(jì)算得到， 也難以描述離子-π相互作用. 根據(jù)相關(guān)實(shí)驗(yàn)數(shù)據(jù)［40］或量子化學(xué)計(jì)算結(jié)果［41］優(yōu)化LJ作用參數(shù)， 是提高力場準(zhǔn)確性的一個(gè)有效途徑， 已被廣泛用于改進(jìn)各類體系［烷烴-水［42，43］、 離子-氨基酸［44，45］、 水-石墨烯［40，46，47］、 氣體-金屬有機(jī)框架材料（MOFs）［48］等］的力場參數(shù). 用此策略發(fā)展離子-石墨烯作用參數(shù)， 可望在保持MD快速、 可模擬較大時(shí)空尺度優(yōu)點(diǎn)的同時(shí)， 準(zhǔn)確描述離子-π作用， 大大提高其模擬石墨烯納米孔道［39，49］的準(zhǔn)確度. 其次， 石墨烯等二維材料與溶液中離子相互作用的熱力學(xué)機(jī)制， 尤其是溶劑對(duì)離子在材料表面吸附的調(diào)控作用目前還缺乏深入研究. 如， 水吸附在石墨烯表面形成密度較高的水層， 對(duì)水分子的氫鍵數(shù)、 居留時(shí)間都會(huì)產(chǎn)生重要影響， 這已得到深入研究［40，46］； 但該水層結(jié)構(gòu)對(duì)離子吸附的影響還少有報(bào)道. 這些問題都限制了二維納米孔道離子輸運(yùn)的研究［7，50］和相關(guān)納米器件的開發(fā)［51］.

本文以石墨烯這一有代表性的（半）金屬二維材料為例， 在Merz團(tuán)隊(duì)［52，53］以及Netz團(tuán)隊(duì)［54，55］開發(fā)的離子力場參數(shù)（分別稱為Merz力場和Netz力場）的基礎(chǔ)上， 優(yōu)化了Li+， Na+， K+， Mg2+， Ca2+和Cl-離子與石墨烯相互作用的LJi-C參數(shù). 優(yōu)化參數(shù)模擬得到的離子-石墨烯吸附自由能（若無特別說明， 只研究溶液中的體系）與量子化學(xué)計(jì)算結(jié)果一致［38］， 可準(zhǔn)確描述（溶液中）石墨烯-離子之間的離子-π相互作用；還將（溶液中）離子在石墨烯表面吸附的平均力勢(shì)［56］（PMFads）分解為石墨烯-離子、 水-離子等相互作用的貢獻(xiàn)（簡稱PMFgra和PMFwat）； 在此基礎(chǔ)上研究了石墨烯表面溶劑結(jié)構(gòu)對(duì)離子吸附的調(diào)控機(jī)制； 進(jìn)一步模擬了實(shí)際濃度下各種鹽溶液/石墨烯界面， 比較了優(yōu)化前后LJi-c參數(shù)對(duì)模擬結(jié)果的影響， 以及這些界面上離子濃度分布與PMFads和離子間相互作用之間的關(guān)系.

1 模型與方法

1.1 模擬系統(tǒng)

模擬了兩種不同類型的體系： 含1個(gè)離子的水/石墨烯界面體系［用于計(jì)算PMFads， 圖1（A）］和鹽溶液/石墨烯界面體系［研究離子的界面行為， 圖1（B）］.

Fig. 1 Simulation system of single ion + graphene in aqueous solution(A) and 1 mol/L salt solution/graphene interface(B)

1.1.1 水/石墨烯界面的單離子 通過傘形采樣［57］模擬計(jì)算（水中）離子在石墨烯表面吸附的PMFads曲線. 如圖1（A）所示， 一個(gè)離子位于C54H18石墨烯薄片（碳、 氫原子均被固定）的中央位置上方， 將其置于含有1351個(gè)水分子的模擬體系（尺寸為3.5 nm×3.5 nm×3.5 nm）. 為了計(jì)算PMFads， 需要進(jìn)行14個(gè)獨(dú)立的模擬； 每個(gè)模擬的離子與石墨烯的垂直（z方向）距離為0.2～1.5 nm（間隔0.1 nm， 各模擬窗口的采樣充分重疊， 表明結(jié)果可靠， 見本文支持信息圖S1）.

在x，y和z方向?qū)﹄x子施加力常數(shù)為1000 kJ·mol-1·nm-2的簡諧勢(shì)以限制離子在特定位置附近運(yùn)動(dòng)：x，y方向簡諧勢(shì)使離子保持在石墨烯中心，z方向簡諧勢(shì)用于傘形采樣. 在正則（NVT）系綜中進(jìn)行20 ns的傘形采樣模擬， 后10 ns模擬軌跡使用加權(quán)直方圖分析方法（WHAM）［58］計(jì)算PMFads曲線. 若無特別說明， PMF計(jì)算均采用碳、 氫原子電荷為0［與圖1（B）體系的原子電荷一致］的C54H18石墨烯片， 體系僅放置單個(gè)離子（無抗衡離子， 這是模擬無限稀釋鹽溶液的通常做法［59］）. 選取PMFads達(dá)到收斂且保持不變時(shí)的勢(shì)能（距離石墨烯大于1.2 nm）作為零勢(shì)能參考面， 取PMFads曲線的最小值［38］（圖S2， 見本文支持信息）為溶液中離子-石墨烯吸附能（Eads）. PMFads的形狀及其最小值（吸附自由能Eads）均隨LJi-C參數(shù)（εi-C）參數(shù)的不同而改變（圖S2）， 通過調(diào)整εi-C參數(shù)， 使得MD模擬得到的Eads與量子化學(xué)方法計(jì)算的水/石墨烯界面的離子吸附能（EQCads， 取自文獻(xiàn)［38］）一致， 從而獲得優(yōu)化后的參數(shù)εion-πi-C（表1）， 參數(shù)εi-C優(yōu)化過程見本文支持信息的S1節(jié).

1.1.2 鹽溶液/石墨烯界面 如圖1（B）所示， 1 mol/L的LiCl， NaCl， KCl， CaCl2或MgCl2鹽溶液（含30個(gè)陽離子與相應(yīng)數(shù)量的Cl-離子， 模擬體系凈電荷為0）與尺寸3.69 nm×3.843 nm的單層石墨烯（540個(gè)碳原子， 電荷設(shè)為0［38，46］，z坐標(biāo)為0）接觸， 置于3.69 nm×3.843 nm×7 nm的模擬盒子. 鹽溶液的厚度約為3.5 nm（含1324個(gè)水分子）， 溶液上方有高度約3.5 nm（模擬盒子高度7 nm）的真空間隙以消除周期性對(duì)模擬結(jié)果的影響. 各鹽溶液/石墨烯界面體系在NVT系綜中模擬20 ns， 采用后10 ns的模擬軌跡進(jìn)行分析.

Table 1 LJi-C parameters and adsorption energy of ion-graphene simulated with Merz and Netz force fields

1.2 模擬設(shè)置

采用LJ 勢(shì)和靜電勢(shì)計(jì)算分子間相互作用； 離子與石墨烯的碳原子之間的LJ 相互作用為V(r)=（其中，εi-C和σi-C是LJi-C參數(shù)）， 通常利用LB 混合規(guī)則得到［εLi-BC=(εi·εC)12，σLBi-C=(σi+σC)/2］.

主要基于量子化學(xué)計(jì)算的（水中）離子在石墨烯表面的吸附自由能Eads［38］優(yōu)化得到離子-碳相互作用強(qiáng)度參數(shù)εion-πi-C（Eads計(jì)算方法與εion-πi-C優(yōu)化方法見1.1.1節(jié)與支持信息的S1節(jié)， 若不特別說明，Eads均指在水中的吸附自由能）； 而σi-C參數(shù)在優(yōu)化前后保持不變（均為σLBi-C）.

分別采用Merz力場［52，53］和Netz力場［54，55］對(duì)離子進(jìn)行模擬； 采用剛性石墨烯模型［石墨烯原子（包括封端的氫原子）的位置均被固定， 這也是石墨烯MD 模擬研究的通常做法［60，61］］， 碳原子力場參數(shù)（εC=0.4899 kJ/mol，σC=0.3214 nm）取自Werder 等［62］的工作； 使用拓展簡單點(diǎn)電荷（SPC/E）水模型［63］. 水、離子以及石墨烯的LJ參數(shù)及原子電荷見表S1（本文支持信息）. 如無特別說明， C54H18的C和H原子電荷均為0； 還計(jì)算了C54H18的限制性擬合靜電勢(shì)（RESP）原子電荷［64，65］， 但原子電荷對(duì)PMF的影響很小，見本文支持信息S2節(jié).

使用GROMACS 5.0.5 軟件［66］進(jìn)行MD 模擬， 使用VMD 軟件［67］進(jìn)行可視化. MD 模擬軌跡采用蛙跳算法進(jìn)行積分， 時(shí)間步長為2 fs. 圖1模擬系統(tǒng)的x，y和z方向都應(yīng)用了周期性邊界條件. 在所有模擬中，分子間LJ相互作用的截?cái)嗑嚯x為1.2 nm，

采用Particle mesh Ewald 方法［68］計(jì)算長程靜電相互作用. 使用SETTLE 算法［69］約束水分子的鍵長. 使用Nose-Hoover方法［70～72］控溫， 耦合溫度為298.15 K. 圖1（A）體系先在NPT系綜進(jìn)行5 ns的預(yù)平衡模擬（使用Parrinello-Rahman 方法［72～74］控壓， 耦合壓力為1×105Pa）， 之后在NVT 系綜中進(jìn)行模擬；圖1（B）的體系使用NVT系綜進(jìn)行模擬.

2 結(jié)果與討論

2.1 離子-石墨烯相互作用參數(shù)的優(yōu)化與離子在水/石墨烯界面的PMFads曲線

研究了Merz和Netz這兩個(gè)團(tuán)隊(duì)分別開發(fā)的Li+， Na+和K+等離子力場參數(shù)（簡稱Merz力場［52，53］、 Netz力場［54，55］， 將未經(jīng)過優(yōu)化的力場參數(shù)簡稱為“優(yōu)化前”力場）. 在這兩套力場參數(shù)的基礎(chǔ)上分別優(yōu)化了離子-石墨烯的作用參數(shù)（簡稱為“優(yōu)化后”力場）： 這些優(yōu)化前后的參數(shù)， 及由其計(jì)算得到的（溶液中）離子在石墨烯表面的吸附自由能Eads見表1.

使用優(yōu)化后力場參數(shù)（εiio-nC-π）計(jì)算得到的Ead（s絕對(duì)值）均比優(yōu)化前力場參數(shù)（εLi-BC）計(jì)算得到的結(jié)果更大， 且與量子化學(xué)計(jì)算結(jié)果（EQCads）［38］一致. 在此基礎(chǔ)上， 測(cè)試了石墨烯片的尺寸大?。▓DS3， 見本文支持信息）、 原子電荷（圖S4， 見本文支持信息）以及抗衡離子（圖S5， 見本文支持信息）等可能影響參數(shù)優(yōu)化的模擬條件（S2節(jié)）， 并計(jì)算相應(yīng)的PMFads曲線（以Na+為例， 使用Merz力場優(yōu)化前、 后的參數(shù)）.

上述模擬條件的改變對(duì)PMFads曲線的影響都很微弱， 對(duì)于C54H18或更大尺寸的石墨烯片， 傘形采樣模擬的Eads與量子化學(xué)計(jì)算的EQCads都基本不變且二者一致［圖S3（A）］， 驗(yàn)證了本文PMFads計(jì)算體系的合理性. 結(jié)果表明優(yōu)化后εion-πi-C參數(shù)可準(zhǔn)確描述溶液中離子-石墨烯之間的離子-π相互作用.

無論用εLBi-C或εion-πi-C參數(shù)模擬得到水/石墨烯界面處離子的PMFads曲線（圖2）， 當(dāng)各離子距離石墨烯z＞ 1.2 nm時(shí)， 它們的PMFads都很接近0， 表明此時(shí)離子-石墨烯相互作用十分微弱； 當(dāng)z＜1.2 nm， PMFads開始劇烈波動(dòng)并產(chǎn)生多個(gè)極大值和極小值， 用εLBi-C或εion-πi-C參數(shù)計(jì)算的PMFads存在明顯差異. 當(dāng)0.3 nm＜z＜1.2 nm，εLBi-C參數(shù)的PMFads基本都大于0（即離子接近石墨烯時(shí)能量上升）， 且PMFads的最小值距離石墨烯較遠(yuǎn)（z＞0.5 nm）， 表明離子不會(huì)吸附在石墨烯表面. 這是因?yàn)槭┑目臻g位阻導(dǎo)致離子水合層變形甚至脫水［40，75］（圖3）， 使PMFads上升的作用處于主導(dǎo)地位（εLBi-C參數(shù)未考慮離子-π作用， 模擬的離子-石墨烯相互作用較弱）. 唯一的例外是Li+離子，εLBi-C參數(shù)的PMFads在石墨烯表面（z＜ 0.5 nm）仍然存在最小值， 但阱深并不深， 與z≈0.75 nm的極小值接近（表1）.

Fig.2 PMFads for Li+(A), Na+(B), K+(C), Ca2+(D), Mg2+(E) and Cl-(F) ions adsorbed at the water/graphene interface

相反，εiio-nC-π參數(shù)的PMFads都比的更低，z＜1.0 nm 時(shí)一般都小于0（僅Li+的個(gè)別距離區(qū)間里，PMFads為很小的正數(shù)）， 這是因?yàn)棣舏io-nC-π參數(shù)可準(zhǔn)確描述離子-π相互作用， 模擬的離子-石墨烯相互作用較強(qiáng)， 使PMFads下降的趨勢(shì)居主導(dǎo)地位. 當(dāng)離子緊密接觸石墨烯時(shí)（z＜0.3 nm）， 無論模擬采用還是εiio-nC-π參數(shù)， 離子的第一和第二水合層都由于石墨烯的空間位阻而脫水（圖3）， 同時(shí)， 石墨烯與離子之間存在LJ勢(shì)能斥力， 因此PMFads曲線均隨z減小而急劇上升.

Fig.3 Relationship between ion-graphene distance and coordination numbers for the first(circle)and second(squares) hydration shells(denoted as Nw1 and Nw2, respectively) for Li+(A), Na+(B),K+(C), Mg2+(D), Ca2+(E) and Cl-(F)

2.2 離子-水和離子-石墨烯相互作用對(duì)PMFads曲線的影響

為了進(jìn)一步研究離子-石墨烯作用和離子-水作用對(duì)PMFads的貢獻(xiàn)， 模擬了以下兩種不同體系：（1） 真空中離子+石墨烯［去除圖1（A）體系的水分子， 其它保持不變］， 離子位于石墨烯正中央的上方；（2） 與圖1（A）模擬體系一致， 但模擬時(shí)不計(jì)算離子-石墨烯作用（設(shè)置εi-C=0 kJ/mol）. 在這兩種體系中，離子分別受石墨烯或水相互作用， 模擬得到的PMF 分別稱之為PMFgra或PMFwat（圖4）， 且PMFads=PMFgra+PMFwat［圖S6（見本文支持信息）， 將PMF分解為模擬各組分貢獻(xiàn)的做法已有研究報(bào)道［76，77］， 表明其準(zhǔn)確性］. PMFgra的形狀比較簡單， 隨著離子接近石墨烯， PMFgra單調(diào)下降，z≈0.3 nm（約等于σi-C，表1）時(shí)達(dá)到極小值（數(shù)值與εi-C成正比， 見圖S7， 見本文支持信息）， 之后， 由于C原子-離子之間的排斥作用急劇上升. 實(shí)際上， PMFgra與離子-石墨烯相互作用的LJ勢(shì)能曲線一致（圖S8， 見本文支持信息）.與之相反， PMFwat的形狀復(fù)雜很多， 當(dāng)離子距離石墨烯比較遠(yuǎn)時(shí)（z＞1.2 nm）， 此時(shí)水可視為各向同性的均勻介質(zhì)， PMFwat是常數(shù). 當(dāng)接近石墨烯時(shí)（z＜1.2 nm）， 一方面水在石墨烯表面形成了兩個(gè)密度較高的水層［32］［圖5（A）］， 這兩層水與離子的作用不同于體相水； 另一方面離子的水合層結(jié)構(gòu)及水合數(shù)均開始變化［33］（圖3、 圖5和圖6）. 這兩方面的作用使得PMFwat有復(fù)雜的變化.

由于離子接近石墨烯時(shí)水合層結(jié)構(gòu)發(fā)生變化（圖3）， PMFwat從z＜1 nm開始總體趨勢(shì)都向上且基本大于0（圖4）； 但陽離子的PMFwat仍存在兩個(gè)極小值（Na+， K+也存在兩個(gè)極小值， 雖較不明顯， 圖4）. 且可見， 陽離子PMFwat極小值的位置與石墨烯表面兩層水的位置［32］存在對(duì)應(yīng)關(guān)系［圖5（A）］： 第一極小值（z=0.3～0.55 nm）位于第一、 第二水層之間； 第二極小值（z=0.7～0.8 nm）位于第二水層之外（遠(yuǎn)離石墨烯）. 另外， 離子第一、 第二水合層開始脫水時(shí)的位置（見本文支持信息表S2， 即水合數(shù)開始變化的位置； 以下分別稱為“第一水合層拐點(diǎn)”和“第二水合層拐點(diǎn)”）分別與PMFwat第一、 第二極小值的位置一致［圖5（B）］. 以上水層位置、 PMFwat極小值位置和水合層拐點(diǎn)之間的關(guān)聯(lián)表明， PMFwat的極小值和極大值可能是離子-水層相互作用和脫水作用兩個(gè)因素共同導(dǎo)致. 進(jìn)一步計(jì)算Li+離子在PMFwat的極小值和極大值位置時(shí)， 離子周圍的水密度分布（圖6， 其它陽離子接近石墨烯時(shí)， 周圍水相對(duì)密度分布變化類似）.

Fig.4 PMFads, PMFgra(contribution of ion-graphene interactions to PMFads) and PMFwat(contribution of ion-water interactions to PMFads) for Li+(A), Na+(B), K+(C), Ca2+(D), Mg2+(E) and Cl-(F) ions adsorbed at the water/graphene interface

Fig.5 Water density distribution above graphene flake(black curve) and the positions of cation PMFwat’s minimum(A), PMFwat’s 1st and 2nd minimum positions agree with corresponding hydration inflection points’ position(B)

Fig.6 PMFwat of Netz force field’s Li+(P1, V1, V2 represent the 1st maximum, the 1st and 2nd minima, respectively)(A), corresponding snapshots of the ion’s 1st hydration shell and the relative density distribution of water(B—D)

為了簡便， 分別用P1， V1 和V2 表示PMFwat的第一極大值和第一、 第二極小值［圖6（A）］. 當(dāng)離子接近第二水層［z≈0.6 nm， 圖5（A）］， 因水層的水密度較高， 離子第二水合層水合數(shù)甚至略有上升［圖3（A）， 由圖6（B）可見， 第二水合層與石墨烯表面水層相交處的水密度增大］， PMFwat開始下降； 當(dāng)離子到達(dá)V2位置［與第二水合層拐點(diǎn)一致， 圖5（B）］， 第二水合層開始脫水（圖3）， 因水合層脫水造成的能量懲罰（若水合層全部脫去， 其能量懲罰即離子的水合能）使得PMFwat上升， 故此時(shí)PMFwat達(dá)到一個(gè)極小值. 隨著離子進(jìn)一步接近石墨烯（z=0.4～0.55 nm）， 石墨烯表面第一水層與離子距離接近［圖6（C）］， 相互作用強(qiáng)烈， 故PMFwat又開始下降并在P1位置產(chǎn)生極大值. 離子在z=0.4～0.6 nm時(shí)處于兩層水中間， 同時(shí)與它們作用［圖5（A）和圖6（D）］， 能量相當(dāng)有利； 而再進(jìn)一步接近石墨烯時(shí)， 離子的第一水合層開始脫水（圖3， 除了Netz力場的Li+）， 使得PMFwat急劇上升（因?yàn)殡x子第一水合層與離子結(jié)合更緊密， 所以PMFwat上升更加劇烈）， 故PMFwat在V1位置達(dá)到另一個(gè)極小值. 進(jìn)一步可見， 各離子的第一水合半徑［離子-（水）氧原子徑向分布函數(shù)的第一極小值位置， 圖S9， 見本文支持信息］與水合層拐點(diǎn)有明顯關(guān)聯(lián)： 水合半徑大的離子， 水合層拐點(diǎn)的位置距離石墨烯也更遠(yuǎn)［圖S10（A）， 見本文支持信息］. 換言之， 水合半徑小的離子需要更接近石墨烯， 其水合層才開始脫水.

離子第二水合層的情況比第一水合層更復(fù)雜， 所以第二水合層拐點(diǎn)與第二水合半徑的關(guān)聯(lián)較弱［圖S10（B）］. 這些結(jié)果說明， 離子接近石墨烯時(shí)， 一方面離子受到石墨烯表面兩個(gè)水層作用， 使PMFwat降低； 另一方面由于離子第一、 二水合層變形甚至脫水使PMFwat升高. 后者一般主導(dǎo)， 使得PMFwat在z＜1 nm 時(shí)呈上升趨勢(shì)， 一般大于0； 而有時(shí)前者也可占主導(dǎo)地位. 因此不同陽離子的PMFwat在z=1 nm以內(nèi)都呈現(xiàn)出一個(gè)極大值和兩個(gè)極小值； 進(jìn)一步接近石墨烯（z＜0.3 nm）時(shí)， PMFwat因離子強(qiáng)烈的脫水作用而急劇上升. 與陽離子相反， 隨著Cl-接近石墨烯， Cl-的PMFwat呈現(xiàn)從水平到單調(diào)上升的趨勢(shì)［圖4（F）］： 在第一和第二水合層拐點(diǎn)之間的區(qū)域（z=0.6～0.7 nm， 表S1）， Cl-的PMFwat基本不變， 因?yàn)榇藭r(shí)Cl-與石墨烯表面水層相互作用的PMFwat下降趨勢(shì)與第二水合層脫水的上升趨勢(shì)相互抵消； 而到達(dá)第一水合層拐點(diǎn)之后（z＜ 0.6 nm）， 脫水作用起主導(dǎo)作用， PMFwat上升.

使用εLBi-C模擬時(shí)， 各陽離子的PMFads與它們的PMFwat都很接近（圖S11， 見本文支持信息）. 在各陽離子PMFads的第一極小值之外， 其PMFgra基本為0， 且PMFgra勢(shì)阱較小無法抵消PMFwat的上升， 所以PMFads與PMFwat很接近. 反之， 使用εion-πi-C模擬時(shí)（圖4）， 由于離子-石墨烯作用強(qiáng)烈， PMFgra不可忽略， 顯著影響了PMFads最小值的大小甚至位置（離子-石墨烯距離）， 因此各離子的PMFads都與PMFwat差距明顯. 進(jìn)一步發(fā)現(xiàn)， Li+， Mg2+和Ca2+離子的PMFads最小值位置與相應(yīng)PMFwat的V1位置一致［圖S12和圖S13（A）（見本文支持信息）Merz力場的Li+的偏差最大， 但不超過0.04 nm］. 而另外3種離子（Na+、 K+和Cl-）的PMFads最小值出現(xiàn)在更靠近石墨烯的位置， 與相應(yīng)PMFgra波谷位置一致［圖S12和圖S13（B）］. 這是因?yàn)椋?Na+， K+和Cl-等離子與水作用比較弱（如， 它們?nèi)軇┗杂赡艿慕^對(duì)值比Li+， Mg2+和Ca2+的?。?2～55］），所以當(dāng)離子接近石墨烯時(shí)， 它們因水合層變形或脫水造成的PMFwat上升比較緩慢， 無法抵消PMFgra的下降趨勢(shì)（考慮離子-π作用后， PMFgra及其變化均不可忽略下降）： 它們的總和（PMFads）也下降， 直到離子進(jìn)一步接近石墨烯， PMFgra開始上升時(shí)為止. 所以Na+， K+和Cl-的PMFads的最小值位置由PMFgra的波谷位置決定（圖S12）. 而Li+， Mg2+和Ca2+則與之相反， 當(dāng)它們水合層變形或脫水時(shí)PMFwat上升迅速， 占主導(dǎo)地位； 故PMFads的最小值位置由PMFwat的V1決定（圖S12）.

2.3 鹽溶液/石墨烯界面上的離子吸附

由于實(shí)驗(yàn)或?qū)嶋H應(yīng)用中與石墨烯接觸的都是一定濃度的鹽溶液. 進(jìn)一步使用εLi-BC及εiio-nC-π參數(shù)模擬了濃度1 mol/L 的鹽溶液/石墨烯界面體系［圖1（B）］， 模擬的電解質(zhì)分別為LiCl， NaCl， KCl， MgCl2和CaCl2］. 使用εLBi-C參數(shù)時(shí)， 各陰陽離子在界面附近的密度分布與體相非常接近（圖7）， 表明離子在石墨烯表面的吸附很微弱. 這是因?yàn)棣臠Bi-C參數(shù)低估了離子-石墨烯相互作用. 相反， 當(dāng)使用εion-πi-C時(shí)， 各離子在石墨烯表面都呈現(xiàn)出明顯的吸附峰（z＞0.5 nm時(shí)迅速衰減至接近體相密度， 圖7）. 這表明強(qiáng)烈的離子-石墨烯相互作用將離子從溶液內(nèi)部吸引到鹽溶液/石墨烯界面上［36，37，78］； 所以能否準(zhǔn)確描述該相互作用對(duì)準(zhǔn)確模擬鹽溶液/石墨烯界面體系至關(guān)重要. 進(jìn)一步發(fā)現(xiàn)， 石墨烯表面各陰陽離子的吸附峰位置與相應(yīng)PMFads最小值的位置十分接近（偏差不超過0.04 nm， 圖S14， 見本文支持信息）. 但是， 各陽離子吸附峰高度與相應(yīng)（單離子）Eads的大小順序并不完全一致； 各鹽溶液的Cl-吸附峰高度也不相同. 此外， 5種鹽溶液-石墨烯界面上， LiCl的Cl-離子吸附峰明顯低于另外4種鹽的情況. 原因可能是， 一價(jià)離子Li+對(duì)Cl-的吸引力不如二價(jià)離子Ca2+和Mg2+， 且在界面上的濃度也低于其它陽離子， 對(duì)Cl-界面吸附的促進(jìn)作用較弱. 這些結(jié)果表明， 在一定濃度的鹽溶液-石墨烯界面體系中， 石墨烯、 水與離子之間的相互作用對(duì)離子濃度分布有重要影響； 而離子-離子之間的相互作用［51］也不可忽略.

Fig.7 Ion’s density profile along 1 mol/L electrolyte solution/graphene interface

3 結(jié) 論

發(fā)展了溶液中Li+， Na+， K+， Mg2+， Ca2+和Cl-離子與石墨烯相互作用的LJi-C力場參數(shù)（優(yōu)化后力場），其模擬得到的離子在（溶液中）石墨烯表面的吸附自由能（PMFads阱深）與量子化學(xué)計(jì)算結(jié)果一致， 表明它們可準(zhǔn)確描述離子與石墨烯之間的離子-π相互作用. 作為對(duì)比， Merz和Netz這兩套力場參數(shù)（優(yōu)化前力場）無法描述離子-π作用， 模擬的離子吸附自由能很接近0， 甚至不吸附. 還將離子在石墨烯表面吸附的PMFads分解為離子-石墨烯（PMFgra）、 離子-水（PMFwat）相互作用的分別貢獻(xiàn). 發(fā)現(xiàn)無論用優(yōu)化前還是優(yōu)化后力場參數(shù)模擬， PMFgra均與離子-石墨烯相互作用勢(shì)能曲線一致. 各陽離子PMFwat在z＜1.0 nm時(shí)都呈現(xiàn)出1個(gè)極大值和2個(gè)極小值（z＞1.0 nm時(shí)基本為0）， 而PMFwat極小值位置、 石墨烯表面兩個(gè)高密度水層的位置、 （離子接近石墨烯時(shí)）水合層開始脫水的位置這三者有明顯關(guān)聯(lián). 進(jìn)一步研究了離子接近石墨烯時(shí)其周圍水密度分布的變化. 這些結(jié)果都說明離子接近石墨烯時(shí)， 其PMFwat由離子水合層變形（或脫水）及其與石墨烯表面水層相互作用這二者共同調(diào)控. 還分別用優(yōu)化前及優(yōu)化后力場參數(shù)模擬了更接近實(shí)際問題的各種1 mol/L鹽溶液/石墨烯界面， 發(fā)現(xiàn)優(yōu)化后力場參數(shù)可準(zhǔn)確描述界面離子吸附現(xiàn)象， 且各陰陽離子的吸附峰位置與對(duì)應(yīng)的PMFads最小值位置基本一致. 作為比較， 優(yōu)化前力場模擬的各陰陽離子在石墨烯表面的吸附非常微弱. 調(diào)節(jié)的是離子與石墨烯的作用參數(shù)， 而不是直接調(diào)節(jié)離子LJ參數(shù)， 所以離子與其它分子（如水）的作用及離子在水溶液中行為均不受影響. 開發(fā)的（溶液中）離子-石墨烯相互作用力場及對(duì)石墨烯表面離子吸附機(jī)制（尤其是其中的溶劑效應(yīng)調(diào)控作用）的探討可為石墨烯或其它二維材料的科學(xué)研究及相關(guān)納米器件開發(fā)提供借鑒和指導(dǎo).

支持信息見http: //www.cjcu.jlu.edu.cn/CN/10.7503/cjcu20230155.

感謝英國曼徹斯特大學(xué)Christopher D. Williams教授關(guān)于溶劑模型和吸附自由能方面的建議.