基于密度泛函理論研究二元排斥Yukawa流體的表面結構性質

1、1460ArticleActaPhys.鄄Chim.Sin.,2006,22(12)：14601465物理化學學報(WuliHuaxueXuebao)December基于密度泛函理論研究二元排斥Yukawa流體的表面結構性質楊摘要振徐志軍楊曉寧*210009)(南京工業(yè)大學化學化工學院,南京基于自由能密度泛函理論(DFT)考察了二元排斥Yukawa(HCRY)流體在不同外場下的密度分布.基于微擾理論,體系的Helmholtz自由能泛函采用硬球排斥部分和長程色散部分貢獻之和,其中Kierlik和Rosinberg的加權密度近似(WDA)被用來計算硬球排斥部分,而色散部分采用平均場理論(MFT)進

2、行描述.為了驗證DFT計硬核直徑和位能參數算結果的合理性,研究中采用巨正則MonteCarlo(GCMC)模擬計算了在不同主體相密度、比的條件下二元HCRY混合流體的密度分布.結果表明,該DFT計算結果與GCMC模擬值吻合良好.關鍵詞：Yukawa流體,中圖分類號：O647密度泛函理論,密度分布,MonteCarlo模擬StructuresofBinaryHard鄄coreRepulsiveYukawaMixturesNearSolidSurfacesfromDFTYANG,ZhenXU,ZhiJunYANG,XiaoNing*210009,P.R.China)(CollegeofChemis

3、tryandChemicalEngineering,NanjingUniversityofTechnology,NanjingAbstractAfreeenergydensityfunctionaltheory(DFT)hasbeenusedtoinvestigatethedensityprofilesofbinaryhardcorerepulsiveYukawa(HCRY)mixturesundertheinfluenceofvariousexternalfields.TherepresentationoftheexcessHelmholtzfreeenergyfunctionalwasba

4、sedontheperturbationtheory,whereaweighteddensityapproximation(WDA)theory,proposedbyKierlikandRosinberg,wasappliedforthehardcorerepulsiveinteractionandameanfieldtheory(MFT)wasusedforthelongrangedispersionone.TotestthevalidityoftheDFT,thegrandcanonicalensembleMonteCarlo(GCMC)simulationswerecarriedoutt

5、oobtainthedensityprofilesofbinaryHCRYmixturesunderdifferentbulkdensities,hardcorediameters,andenergyparameterratios.TheDFTcalculationsareinquitegoodagreementwiththeGCMCsimulations.Keywords：Yukawafluid,Densityfunctionaltheory,Densityprofile,MonteCarlosimulation最近,作為一種分子模型的Yukawa勢能函數受到廣泛的關注,此模型已經被成功用于

6、模擬簡單流體1、膠體23、等離子體45、C60分子67和球型蛋白質分子810.文獻中已有許多有關色散作用為Yukawa熱力學以及結構形式的流體的相平衡、界面張力、性質的研究報道.付東等11基于熱力學二階微擾理論,結合密度泛函理論(DFT)、DuhMierYTeran狀態(tài)方程和直接相關函數的解析解,建立了適用于均Received:June2,2006;Revised:July7,2006.國家自然科學基金(20476044)資助項目*勻與非均勻流體的狀態(tài)方程,研究了三種不同范圍參數姿的純Yukawa流體的汽液相平衡和界面張力.Tang12利用一階平均球近似(MSA)和Rosenfeld的微擾方法

7、對純Yukawa流體在硬墻和吸引墻表面上的密度分布進行了研究,得到了令人滿意的結果.積分方程理Cochran和Chiew13分別采用微擾理論、論和MonteCarlo模擬方法對純Yukawa流體的熱力學性質和結構性質進行了研究.Zhou和Jamnik14Correspondent,Email：Yangxia;Tel:02583587185.鬁EditorialofficeofActaPhysicoChimicaSinicaNo.12楊曉寧等：基于密度泛函理論研究二元排斥Yukawa流體的表面結構性質1461根據二階直接相關函數中的參考項與微擾項對主體相密度依賴程度不同,分別采用三階和二階泛函微

8、擾展開近似(functionalperturbationexpansionapproximation,FPEA),并通過該微擾DFT研究了純Yukawa流體在三種不同外場存在時的結構性質.Yu等15采用修正的基本度量理論(modifiedfundamentalmeasuretheory,MFMT)和Rosenfeld的微擾方法,并通過一階MSA得到主體相的直接相關函數,以此研究了純Yukawa流體在狹縫孔中的結構性質以及超額吸附,并與相應的模擬值進行了對比,結果表明理論計算結果是令人滿意的.然而目前文獻中對于更為復雜的二元Yukawa混合流體的研究相對較少.You等16通過簡單地在平均場理論(

9、MFT)中引入一個修正參數c來描述長程色散作用對自由能的貢獻,同時采用MFMT表達硬球作用對自由能的貢獻并由此建立了MFMTcMF理論,考察了二元Yukawa混合流體在狹縫孔中的結構性質.隨后,Patra和Smith17基于均勻體系得到的自由能泛函形式用于描述相應的非均勻體系,發(fā)展了一種具有熱力學一致性的DFT方法,該方法計算得到的結果和參考文獻16中報道的模擬值吻合良好.以上提及的二元Yuakwa流體的研究報道都是色散作用為吸引形式的Yukawa(hardcoreattractiveYukawa,HCAY)流體.然而據我們所知,對于色散作用為排斥形式的Yukawa(hardcorerepul

10、siveYukawa,HCRY)流體的研究均是關于單組分的,有關二元HCRY混合流體的研究至今未見文獻報道.密度泛函理論作為一種倍受關注的分子理論計算方法,因其具有良好的理論基礎和計算簡單的優(yōu)勢,已經成為研究非均勻體系的重要方法.DFT方法處理的關鍵問題在于如何來確定過量Helmholtz自由能的泛函形式.基于微擾理論,通常將其分為硬球排斥項和長程項部分.對于硬球排斥項部分一般有局部密度近似(localdensityapproximation,LDA)和加權密度近似(weighteddensityapproximation,WDA)兩種處理方法.對于混合物的WDA方法中,最著名的是Rosenf

11、eld18提出的基本度量理論(fundamentalmeasuretheory,FMT).一方面Roth19和Yu2021等基于BMCSL方程對FMT進行了改進,克服了FMT對壁面接觸點處密度的預測值偏高這一缺點.另一方面Kierlik和Rosinberg2223對FMT中的權函數重新進行了設計,用標量權函數代替了FMT中表示表面密度不均勻的矢量權函數.同時還證明了24FMT的權函數與他們所設計的更為簡單的權函數形式在本質上是等價的,而由于這種權函數形式與密度本身沒有關系,使得計算得到了很大的簡化而被普遍采用.本文的主要目的是采用Kierlik和Rosinberg的DFT來研究非均勻二元HCR

12、Y混合流體在固體表面的結構性質,以及通過巨正則MonteCarlo(grandcanonicalMonteCarlo,GCMC)模擬方法驗證了在不同條件下該理論計算結果的合理性.此外,為了更好地了解Yukawa流體分子間的長程色散作用對其結構性質的影響,計算比較了相同條件下二元HCAY混合流體與二元HCRY混合流體的密度分布.11.1理模論型HCRY混合流體分子間相互作用勢表示為,u,r<ijij(r)=,rij(1)式中ij=(i+j)/2,ij=(ij)1/2,i,j=1,2;i是組分i的硬核直徑,i為兩分子間作用的位能參數,r是兩分子間的核心間距,是Yukawa勢能作用的屏蔽長度參

13、數,=1.8.=1/kBT,kB為Boltzmann因子,T為絕對溫度.本文采用對比溫度,T*=kBT/1=4.同時,HCRY流體分子與固體表面的相互作用勢可用下式來表示：Vexti(z)=wexp(-(z-i/2)/i),zi,otherwise(2)式中z為流體分子與墻面的垂直距離.當w=0時,式(2)代表硬球作用的表面;w<0時,式(2)代表吸引作用的表面;當w>0時,式(2)代表排斥作用的表面.1.2密度泛函理論流體密度泛函理論是以熱力學函數中的巨勢作為起點,巨勢可以表示成局部密度i(r)的唯一泛函形式.基于微擾理論,HCRY流體分子間的相互作用可以分成硬球排斥項和長程色散

14、項,uij(r)=uhsij(r)+udisij(r).那么,過量的Helmholtz自由能泛函也被分成硬球排斥和色散部分,Fex=Fexhs+Fexdis.據我們所知,對于二元HCRY混合流體的二階直接相關函數并沒有解析的表達式,因而采用MFT對色散部分的自由能泛函進行求解.對于一個給定的體系,巨勢可以表示為,1462ActaPhys.鄄Chim.Sin.(WuliHuaxueXuebao),2006Vol.222i(r)=kBTi=1ri(r)ln(i(r)3i-1+Fexhs2di(r)+i=122dri(r)Vexti(r)-i+drdr(r)(r)udisi=1j=1ijij(|r-

15、r|)(3)式中i為deBroglie熱波長,i為組分i的化學勢.過量的Helmholtz自由能泛函的硬球排斥項Fexhs在Kierlik和Rosinberg2223的WDA方法中表示為,Fexhsi(r)=kBTdrn(r)(4)式中是均勻硬球流體的Helmholtz自由能密度,n(r)是空間位置r處的權重密度,可以表示為,n2(r)=dri(r)i()(|r-r|)(5)r分i=1式中r、別表示不同位置處的位置矢量.按照定標粒子理論(scaledfieldparticletheory,SPT),能夠表示成如下形式2223：=-n0ln(1-n3)+3+3(6)式中n0、n1、n2和n3為S

16、PT中的對比變量：n2=i=1iRi()(7)以及：Ri(0)=1,Ri(1)=Ri,Ri(2)=4R2i,Ri(3)=R3i(8)式中Ri為組分i的硬核半徑.四個權重函數i()(r)簡單地與Heaviside階梯函數(r)的逐次導數(、和)相關2223：i(3)(r)=(Ri-r)(91)i(2)(r)=(Ri-r)(92)i(1)(r)=(Ri-r)(93)i(0)(r)=-(Ri-r)+(Ri-r)(94)通過對體系的巨勢最小化來求得平衡時的密度分布:i(r)/i(r)=0(10)2巨正則MonteCarlo模擬為了驗證KRDFT理論計算結果的合理性,我們進行了一系列GCMC模擬(i、V

17、、T不變)計算.模擬盒子為一個邊長101的正方體.Yukawa勢能函數的截斷半徑取51,超過此距離,Yukawa流體間相互作用可以忽略不計.HCRY混合流體受限于垂直z軸方向的兩平行墻之間.兩墻以原點為中心,墻面間距為101,此距離也足夠使得兩墻之間不發(fā)生勢能的重疊.同時僅僅在x和y軸方向上采用周期性邊界條件.對于二元混合體系,標準的GCMC模擬方法中每一步由四種基本嘗試組成25:(1)隨機選取一個粒子,嘗試隨機移動到另一個空間位置;(2)隨機選取一個空間位置,嘗試插入一個粒子;(3)隨機選取一個粒子,嘗試刪除該粒子;(4)隨機選取一個粒子,嘗試將該粒子變換成另一種粒子.其中最大位移的調整以使

18、得嘗試(1)的接受概率為50%左右為準.以上四種嘗試在模擬中采用2111的比例進行,每一種嘗試中兩種粒子等概率被選取.在每個條件下,GCMC模擬總共運行1.2×108步,前2.0×107步用來平衡.平衡后,每隔500步保存一次,并進行統(tǒng)計.此外,為了使得GCMC和DFT兩者計算條件保持一致,采用如下方法26得到相應主體相密度(bi,i=1、2).基于相同的輸入條件(i、V和T),對主體相流體亦進行了GCMC模擬計算.在x、y和z軸方向上均采用周期性邊界條件.GCMC中的移動、插入、刪除和變換四種嘗試在模擬中采用111:1的比例進行,每一種嘗試中兩種粒子等概率被選取.在每個條

19、件下,GCMC模擬總共運行4.0×107步,前1.0×107步用來平衡.平衡后,每隔200步保存一次,并進行統(tǒng)計,同時將得到的主體相密度(bi,i=1、2)作為DFT的輸入條件.3.1結果與討論二元HCRY混合流體首先,我們采用GCMC研究了不同條件下,二元HCRY混合流體在固體表面的密度分布.此外,通過對相應條件下的主體相流體模擬計算,得到該條件下的主體相密度(bi,i=1、2),然后將該主體相密度(bi,i=1、2)作為輸入條件進行DFT計算.為了清晰起見,不同條件的體系參數及主體相密度匯總在表1中.圖1與圖2分別顯示了在條件I和條件II下,二元HCRY混合流體在固體表

20、面的密度分布.從圖中我們可以看到,當體系的主體相密度較低(b=b131+b231=0.197)時,如圖1所示,流體密度在固體表面附近處未出現密度的振蕩.當流體接近固體表No.12楊曉寧等：基于密度泛函理論研究二元排斥Yukawa流體的表面結構性質1463表1不同條件的體系參數Table1ThesystemparametersunderdifferentconditionsSystem2/12/1b131b231wBII1.00.80.1180.335-1.001.0III2.00.60.5410.120-1.01.0面時,流體密度分布開始受到固體表面的作用,圖1(a)中的表面對流體分子具有吸引

21、作用,使得流體在表面附近的密度逐漸地升高,出現了局部密度大于體相密度的富集區(qū),呈現出正吸附狀態(tài);而表面為排斥作用時,如圖1(c)所示,流體在表面附近的局部密度非但沒有升高,相反較體相密度有細微的下降,呈現出負吸附狀態(tài).這是因為表面的排斥作用使得流體分子遠離表面,大大削弱了流體在表面的堆積程度所致.隨著體系的主體相密度(b=b131+b231=0.453)逐漸升高,如圖2所示,流體分子間的體積排斥效應也逐漸增大,使得流體分子更易靠近固體表面.因而流體在表面的堆積程度也逐漸增大以及在圖1在條件I下二元HCRY混合流體在固體表面的密度分布Fig.1Thedensityprofilesofbinary

22、HCRYmixturesnearsolidsurfacesforthesystemIThesymbolsandlinesarefortheGCMCandDFTresults,respectively.w/(kBT):(a)-1.0,(b)0,(c)1.0表面的接觸點處達到了更高的密度峰值.尤其在圖2(c)中,流體分子間體積排斥效應的增大使得流體呈現出與圖1(c)中正好相反的正吸附狀態(tài).在條件I和II下,DFT的計算結果與模擬值吻合良好,只是條件II中兩者在表面的附近存在著細微的偏差.圖3中對比了在條件III下二元HCRY混合流體密度分布理論值與模擬值.從圖中我們可以看到,理論值與模擬值在接觸點

23、附近存在著一定的偏差,表面為吸引作用時,兩者偏差最大.但就整體而言,DFT的計算結果較好地吻合了GCMC的模擬結果.在體系III中,主體相密度(b=b131+b231=0.661)處在一個較高的位置,導致流體分子間表現出強烈的體積排斥效應.流體密度在表面附近處陡然地升高,此時體系的密度分布曲線呈現出明顯的振蕩衰減形式.對于組分1而言,流體在表面的接觸點(z/1=0.5和z/1=1.6)處各有一個峰值.兩個波峰的距離相差1.11,要大于其本身的硬核直徑1,這與文獻中報道27純HCRY流體的密度分布曲線有更長的波圖2在條件II下二元HCRY混合流體在固體表面的密度分布Fig.2Thedensity

24、profilesofabinaryHCRYmixturenearsolidsurfacesforthesystemIIThesymbolsandlinesarefortheGCMCandDFTresults,respectively.w/(kBT):(a)-1.0,(b)0,(c)1.01464ActaPhys.鄄Chim.Sin.(WuliHuaxueXuebao),2006Vol.22圖3在條件III下二元HCRY混合流體在固體表面的密度分布Fig.3ThedensityprofilesofabinaryHCRYmixturenearsolidsurfacesforthesystemIII

25、ThesymbolsandlinesarefortheGCMCandDFTresults,respectively.w/(kBT):(a)-1.0,(b)0,(c)1.0動周期是一致的.當流體離表面的距離增大到一定程度時,密度分布曲線又變成一條直線,回歸到均勻體系的形式.與圖1和圖2所示的一樣,表面為吸引作用時,密度分布有著最為明顯的振蕩;而對于具有硬球和排斥作用的表面,振蕩的幅度依次減弱.3.2二元HCAY混合流體為了更好地比較二元HCAY與二元HCRY混合流體結構性質的區(qū)別,采用DFT以及GCMC計算了在條件II下,二元HCAY混合流體在固體表面的密度分布.如圖4所示,相比于HCRY流體,

26、HCAY流體不易在表面附近形成堆積,同時在表面的接觸點處的密度峰值也較低.這種區(qū)別在Zhou和Jamnik14研究純HCRY流體的結構性質中也有同樣的結論.尤其當表面為排斥作用時,如圖4(c)所示,HCAY流體出現了明顯負吸附現象,而相同條件HCRY流體則表現出正吸附狀態(tài).兩種流體的密度分布曲線表現出截然不同的形式,主要是因為流圖4在條件II下二元HCAY混合流體在固體表面的密度分布Fig.4ThedensityprofilesofabinaryHCAYmixturenearsolidsurfacesforthesystemIIThesymbolsandlinesarefortheGCMCan

27、dDFTresults,respectively.w/(kBT):(a)-1.0,(b)0,(c)1.0體分子間的不同相互作用力造成的.HCAY流體分子間表現出的相互吸引作用使得流體分子遠離墻面,此外,表面的排斥作用同樣使得流體分子遠離表面,這種雙重作用使得流體在表面呈現出更為明顯的負吸附狀態(tài).從圖4中還可以看到,理論值與模擬值除在緊貼表面的附近處存在著較小的偏差外,整體上兩者吻合良好.說明在條件II下,該DFT對二元HCAY混合流體的計算結果是令人滿意的.然而,關于Kierlik和Rosinberg能否合理地描述二元HCAY流體的結構性質還有待進一步的研究.3.3組分的過量吸附量為了定量地描

28、述流體在固體表面的吸附情況,結合主體相密度和平衡密度分布,組分i在固體表面的過量吸附量exi可由下式計算得到：exi=不同條件下乙肄i/2(i(z)-bi)dz(11),由上式積分計算得到的各組分的表面No.12楊曉寧等：基于密度泛函理論研究二元排斥Yukawa流體的表面結構性質1465表2不同條件下組分的過量吸附量Table2ThesurfaceexcessforeachcomponentunderdifferentconditionsSystemI(HCRY)GCMC0.04080.0138-0.00940.05090.0149-0.0131II(HCRY)DFT0.04950.03020

29、.00780.12710.06470.0089GCMC0.04460.02860.00880.11520.06890.0103II(HCAY)DFT0.03180.0016-0.02880.09730.0106-0.1754GCMC0.03110.0020-0.02780.10070.0142-0.1668III(HCRY)DFT0.17210.14650.11270.05810.03230.0171GCMC0.15060.13020.10080.04520.02550.0192wB過量吸附量匯總在表2中.通過對比表2中DFT的計算值與GCMC的模擬值,表明DFT的計算結果較為準確地預測了三個

30、條件下各組分的過量吸附量.流體分子間的體積排斥效應以及固體表面的吸引作用均使得流體分子靠近表面,從而增大了流體在固體表面的吸附程度.相反地,固體表面的排斥作用卻抑制了流體在固體表面的吸附.正如表2中所示,隨著體系密度的升高,流體分子間的體積排斥效應增大,流體在固體表面的過量吸附量也逐漸地增大.當表面為吸引作用時,此時各組分的過量吸附量最大,而表面為排斥作用時,則各組分的過量吸附量最小.此外,在條件II下,HCRY流體較HCAY流體更易被吸附,有著更大的過量吸附量.4結論本文采用Kierlik和Rosinberg的DFT結合平均場理論計算了在固體表面二元HCRY混合流體的密度分布及各組分的過量吸

31、附量.流體的密度分布及其在固體表面的吸附均是由流體分子間以及流體分子與固體表面的相互作用所共同決定的.隨著體系密度逐漸升高,流體分子間的體積排斥效應增大,使得流體的密度分布曲線呈現出更為強烈的振蕩衰減形式,同時表面的過量吸附量也逐漸增大.相比于HCAY流體,HCRY流體分子間的排斥作用使得流體分子更易在表面形成堆積,從而有著更為明顯的密度振蕩和更大的過量吸附量.此外,通過與GCMC模擬結果比較表明,當體系密度較低時,理論結果很好地吻合了模擬值;隨著體系密度升高,在表面附近處兩者呈現出一定的偏差,但就整體而言,理論值與模擬結果仍然吻合良好.因此,Kierlik和Rosinberg的DFT方法能夠

32、成功地描述二元HCRY混合流體的表面結構性質.References1Tang,Y.P.;Tong,Z.F.;Lu,B.C.Y.FluidPhaseEquilibria,1997,134:212Davoudi,B.;Kohandel,M.;Mohammadi,M.;Tanatar,B.Phys.Rev.E,2000,62:69773Fu,D.;Li,Y.G.;Wu,J.Z.Phys.Rev.E,2003,68:0114034Totsuji,H.;Kishimoto,T.;Totsuji,C.Phys.Rev.Lett.,1997,78:31135VanHorn,H.M.Science,1991,252:3846Lomba,E.;Almarza,N.G.J.Chem.Phys.,1994,100:83677Germain,P.;Amokrane,S.Phys.Rev.E,2002,65:0311098Lin,Y.Z.;Li,Y.G.;Lu,J.F