廣義等溫等壓系綜-分子動力學(xué)模擬全原子水的氣液共存形貌?

尹靈康 徐順Seongm in JeongYongseok Jho 王健君 周昕?

1)(中國科學(xué)院大學(xué)物理科學(xué)學(xué)院,北京 100049)

2)(中國科學(xué)院計算機網(wǎng)絡(luò)信息中心,北京 100190)

3)(Center for Soft and Living M atter,Institute for Basic Science(IBS),U lsan 44919,Repub lic of Korea)

4)(中國科學(xué)院化學(xué)研究所,北京 100190)

廣義等溫等壓系綜-分子動力學(xué)模擬全原子水的氣液共存形貌?

尹靈康1)徐順2)Seongm in Jeong3)Yongseok Jho3)王健君4)周昕1)?

1)(中國科學(xué)院大學(xué)物理科學(xué)學(xué)院,北京 100049)

2)(中國科學(xué)院計算機網(wǎng)絡(luò)信息中心,北京 100190)

3)(Center for Soft and Living M atter,Institute for Basic Science(IBS),U lsan 44919,Repub lic of Korea)

4)(中國科學(xué)院化學(xué)研究所,北京 100190)

(2017年3月29日收到;2017年5月4日收到修改稿)

探測相變過程中瞬時共存相的形貌等特征對理解其微觀機制十分重要.本文應(yīng)用廣義等溫等壓系綜-分子動力學(xué)模擬方法,研究全原子水模型的氣液兩相平衡及相變的中間過程.研究發(fā)現(xiàn),此廣義系綜方法能夠通過持續(xù)降溫,連續(xù)地歷經(jīng)從氣態(tài)、氣液共存到液態(tài)的整個相變過程,通過持續(xù)升溫歷經(jīng)其相反過程,而不會發(fā)生標(biāo)準(zhǔn)正則系綜中的過飽和熱滯現(xiàn)象.該方法不需要使用副本交換等增強抽樣方法,因此可以用于較大體系的研究,多個獨立的模擬即能獲得整個氣態(tài)液態(tài)區(qū)的平衡性質(zhì)及共存相特征.本文還提出了計算氣液共存界面面積的新方法,給出了水的氣液共存界面形狀隨溫度、壓強的變化規(guī)律.結(jié)果表明,低壓時水的氣液共存界面因其較大的表面張力而接近球面,符合經(jīng)典成核理論的描述,但隨著壓強升高接近其臨界壓強時,氣態(tài)和液態(tài)的差別減小,界面的表面張力變小,界面形狀變?yōu)闊o規(guī)則的枝杈結(jié)構(gòu),表現(xiàn)為二階相變特征.

廣義正則系綜,氣液相變,氣液共存,分子動力學(xué)

1 引 言

氣液相變是自然界中非常普遍的現(xiàn)象,在人類的生活生產(chǎn)中扮演著重要角色,然而人們對相變的微觀機制并沒有完全理解.描述相變過程的一般理論是經(jīng)典成核理論(CNT)[1-3].CNT認(rèn)為相變的發(fā)生受成核過程支配,即形成一個具有一定大小的、球形的新相的核,稱為臨界核,然后這個臨界核迅速長大,直至整個系統(tǒng)完全變成新相.但是CNT對臨界核的形狀、結(jié)構(gòu)做了過于簡單的近似,不能解釋很多實驗現(xiàn)象.近來,大量工作著眼于研究相變過程中的微觀細(xì)節(jié)[4-8],人們發(fā)現(xiàn)用多步成核理論[1,9-12]來描述許多相變成核過程可能更加精確.這些結(jié)果表明,想要清晰地理解相變過程,形成核以及核生長變大的微觀細(xì)節(jié)是十分重要的.

人們發(fā)展了許多方法來研究相變成核過程[13-17]以及相共存態(tài)的微觀結(jié)構(gòu)[18-23]. Zahn[24]采用路徑抽樣技術(shù),通過分子動力學(xué)模擬研究簡單模型的水的沸騰現(xiàn)象,給出了液態(tài)體系中相變開始階段氣相核的形成過程.Panagiotopoulos[25]于1987年提出了Gibbs系綜Monte Carlo(GEMC)方法,可以直接獲得共存區(qū)的相平衡.M ou?ka和Nezbeda[26]以及Trejos等[27]應(yīng)用該GEMC方法分別研究了BK模型水的氣液平衡過程和受限量子液體的氣液共存態(tài).Cho等[28-31]對勢能函數(shù)做出修改,使用廣義副本交換方法,能夠抽樣到穩(wěn)定存在的共存區(qū)的構(gòu)象,研究了Lennard-Jones流體和全原子水等多種模型的兩相共存形貌特征.Xu等[32]引入廣義正則系綜(GCE),結(jié)合副本交換技術(shù),使用Monte Carlo (MC)方法模擬了格點Potts模型的鐵磁順磁相變共存區(qū).Jeong等[33]應(yīng)用這個GCE副本交換模擬,給出了Lennard-Jones模型和粗?；痬W水模型的氣液共存相的微觀形貌隨溫度和壓強的變化.

這些相變共存態(tài)的模擬研究通常使用副本交換等增強抽樣技術(shù),因此一般集中在較小或較簡單的系統(tǒng).本文擴展GCE的MC模擬到分子動力學(xué)模擬(MD),研究大尺寸全原子水體系的氣液轉(zhuǎn)變及其共存相的微觀結(jié)構(gòu).相較于MC模擬,MD方法有許多成熟的軟件包可供使用,因此在研究全原子水模型等比較復(fù)雜的系統(tǒng)中比較方便.我們在MD通用程序包(LAMMPS,GROMACS等)中實現(xiàn)了廣義等溫-等壓系綜MD方法(gNPT-MD).

我們模擬了4000個TIP4P/2005水分子在不同壓強下的氣液相變,發(fā)現(xiàn)至少在壓強不太低的情況下,等溫-等壓(NPT)系綜下水的氣液相變的過飽和熱滯現(xiàn)象在廣義等溫-等壓(gNPT)系綜中被很好地消除,系統(tǒng)能連續(xù)地在氣液兩相之間轉(zhuǎn)變.從而無需使用副本交換等技術(shù),通過持續(xù)升溫或降溫的gNPT-MD模擬,我們就可以得到系統(tǒng)在整個相變區(qū)的微觀構(gòu)象,包括純液態(tài)、過熱液態(tài)、氣液共存、過冷氣態(tài)以及純氣態(tài)等,并準(zhǔn)確地給出了氣態(tài)液態(tài)間的相平衡條件.為了定量分析共存相的微觀特征,我們嚴(yán)格定義液滴表面分子可達面積為液氣共存界面面積,提出了一個較精確地計算氣液兩相界面面積的方法,得到了共存區(qū)內(nèi)氣液界面的形貌隨溫度、壓強的變化規(guī)律.對比不同壓強下共存區(qū)的構(gòu)象,我們發(fā)現(xiàn)低壓時的共存界面更接近于球面,基本符合經(jīng)典成核理論,而高壓時界面形狀復(fù)雜,存在許多枝杈結(jié)構(gòu),表現(xiàn)為二階相變的特征.計算得到的表面張力表明,高壓時表面張力較小,不足以克服界面的形狀漲落.

本文的組織如下：第2部分先簡單地介紹了GCE的原理,gNPT-MD在模擬軟件包中的實現(xiàn),以及gNPT系綜和NPT系綜的關(guān)系,然后給出了Voronoi單胞分析方法,并詳細(xì)說明我們改進的計算共存界面面積的方法,隨后給出模擬的細(xì)節(jié);第3部分,給出了用gNPT-MD方法來研究全原子水模型的氣液相變及共存區(qū)界面形貌、表面張力等主要結(jié)果;第4部分給出結(jié)論.

2 研究方法

2.1 GCE與gNPT系綜

相變共存區(qū)在常規(guī)系綜(如正則系綜NVT, NPT等)下不穩(wěn)定,其微觀構(gòu)象不能被有效探測.廣義系綜能夠通過使用能量依賴的熱源溫度(或等價地修改勢能面)而有效地抽樣這些共存構(gòu)象.

NVT系綜的構(gòu)象分布函數(shù)可以寫為

這里β0是倒溫度,E(r)是系統(tǒng)勢能,r是構(gòu)象坐標(biāo)的簡單記號.在廣義系綜如GCE中,勢能E(r)被有效勢能Ug(E)代替：

其中β0,E0和α為常數(shù)參數(shù).改變β0或E0可得到不同的GCE.α＞0對應(yīng)GCE,α=0對應(yīng)NVT系綜.

在分子動力學(xué)模擬中實現(xiàn)GCE系綜非常簡單,只需在標(biāo)準(zhǔn)的NVT系綜模擬中將熱源溫度替換為能量依賴的瞬時溫度 T(E)=1/βg(E)[32],這里GCE的倒溫度

無需修改模擬程序的其他部分.這里E為當(dāng)前構(gòu)象的勢能,因此GCE中熱源溫度在模擬的每一步均被重置.容易證明,這樣可以得到GCE下的構(gòu)象空間分布函數(shù),P(r)∝ e-Ug(E(r)).值得注意的是,這樣實現(xiàn)的GCE模擬僅保證構(gòu)象空間內(nèi)等式(2)給出的這個有效勢能的分布,其在動量空間的分布不是任何溫度下的Maxwell-Boltzmann分布.

在gNPT系綜中,等式(2)中的勢能E用焓H=E+PV代替,有效勢能可以寫為

這里P和V分別是外部壓強和系統(tǒng)體積,H0是常數(shù).同樣可以有g(shù)NPT系綜中焓依賴的熱源倒溫度：

在任意標(biāo)準(zhǔn)分子動力學(xué)軟件包中,實現(xiàn)gNPT系綜與實現(xiàn)GCE類似,即通過等式(5)計算當(dāng)前構(gòu)象的H對應(yīng)的熱浴的瞬時溫度.值得注意的是,在計算系統(tǒng)當(dāng)前構(gòu)象壓強以耦合目標(biāo)壓強時,系統(tǒng)壓強中理想氣體部分的貢獻由當(dāng)前重置后的熱浴溫度直接給出,不使用系統(tǒng)當(dāng)前的動能來計算,而系統(tǒng)瞬時壓強中的維里部分無需任何改變.

2.2 gNPT系綜與NPT系綜的關(guān)系

我們以gNPT與NPT系綜為例說明廣義系綜與正則系綜的關(guān)系,GCE和NVT系綜的關(guān)系基本類似.任意系統(tǒng)的構(gòu)象空間性質(zhì)由其態(tài)密度eS(H),或者內(nèi)稟倒溫度函數(shù)

確定,這里S(H)是熵函數(shù).對gNPT系綜,熱源倒溫度βgNPT(H)和內(nèi)稟倒溫度曲線βs(H)的交點給出H空間分布的極值點.如果此處βgNPT(H)的斜率比βs(H)的切線的斜率大,則為分布的極大值點,即該gNPT系綜主要訪問這個附近.對較大的α,這個極大值點惟一,近似等于此gNPT系綜下H的平均值.

圖1給出了典型的相變系統(tǒng)的βs(H)曲線.對NPT系綜,熱源倒溫度與H無關(guān),對應(yīng)斜率為0的水平直線,而gNPT系綜的熱源倒溫度為斜率為α的直線.在相共存溫度區(qū)間,每個NPT系綜的倒溫度線與 βs(H)有三個交點,如圖中A,B和C點.但僅純相或過飽和相(A和C)是概率分布的極大值點,能在此NPT系綜中被充分訪問,而共存相B是分布的極小點,僅有非常小的訪問概率,歸因于B點處βs(H)曲線的正斜率.在NPT系綜中,在降溫和升溫過程中,系統(tǒng)分別從氣態(tài)到液態(tài)以及從液態(tài)到氣態(tài)的相變的βs(H)-H曲線不重合,而是存在熱滯現(xiàn)象,如圖中帶箭頭的紅色虛線所示.因此通常的NPT模擬很難準(zhǔn)確地給出兩相共存溫度,難以有效地探測兩相共存的焓區(qū)(內(nèi)稟倒溫度的正斜率區(qū))的微觀構(gòu)象.而這些共存構(gòu)象對應(yīng)不同溫度的NPT系綜中相變發(fā)生的過渡態(tài),對理解相變微觀機制十分關(guān)鍵.在gNPT系綜中,當(dāng)熱源倒溫度函數(shù)的斜率α大于系統(tǒng)內(nèi)稟倒溫度曲線在共存區(qū)域的斜率時,每個gNPT系綜的倒溫度線與βs(H)只有一個交點,因此共存構(gòu)象能被充分抽樣.

可以給出系統(tǒng)內(nèi)稟倒溫度曲線 βs(H)上的一點(,).改變β0(或者H0)的值,gNPT系綜可以充分抽樣系統(tǒng)任意焓區(qū)的微觀構(gòu)象,給出整個βs(H)曲線,并能夠積分得到熵函數(shù)S(H).

圖1 (網(wǎng)刊彩色)典型的相變過程的倒溫度曲線βs(H)與自由能曲線 最上方的藍色曲線是相變過程的倒溫度曲線,紅色虛線是NPT系綜給出的相變曲線,紅色直線為NPT系綜和gNPT的熱源倒溫度線,下方三條藍色曲線分別給出了T1,Tc,T2三個溫度下的自由能曲線(β=1/kB T)Fig.1.(color on line)Typical cu rves of the recip rocal of temperature of phase transitionβs(H)and cu rves of free energy:The b lue cu rve on the top rep resent βs(H)of typical phase transition,the red dashed lines correspond to phase transition of NPT ensem b le,the red fu ll lines correspond to the recip rocal of temperature of NPT ensem b le and gNPT,and the other b lue curves show the free energy at T1,Tc,T2(β=1/kB T).

由此我們?nèi)菀椎玫饺我鉁囟萒下NPT系綜的H空間分布,PNPT(H)∝e-H/T+S(H),以及自由能,

這里S(H)由變化參數(shù)β0而得到的一系列g(shù)NPT系綜下的平均焓及對應(yīng)的給出的βs(H)曲線積分得到,

S(H)依賴于系統(tǒng)外部壓強P.對每個壓強下的曲線 βs(H),我們可以確定共存溫度 Tc,使得βc=1/Tc與 βs(H)所包圍的上下兩部分面積相等,即圖1中兩個陰影部分面積相等.此時氣液兩相的自由能相等,自由能曲線的能壘由表面能提供.一般地,對任意溫度和壓強,自由能能壘為,

這里cg和cl分別是能壘H?處氣態(tài)和液態(tài)分子的比例,可由H?=cgHg+clHl以及cg+cl=1給出;Hg和Hl分別是此溫度、壓強下氣態(tài)和液態(tài)的平均焓值;γ(T,P)是界面張力,A是界面面積.我們可以由此計算出界面的表面能及表面張力[34,35].

2.3 Voronoi單胞分析方法

Voronoi單胞[36-38]是一種重要的分析工具,可以解決一些隨機介質(zhì)中的結(jié)構(gòu)分析問題,如玻璃、多胞固體、蛋白質(zhì)等.把空間中任意一點與其他點相連,所有連線的垂直平分面可以包圍成一個封閉的多面體,這個多面體稱為Voronoi單胞. Voronoi單胞可以把一個封閉空間劃分成許多充滿空間、但不交疊的凸多面體.在本文中,我們通過Voronoi單胞分析方法來分析模擬結(jié)果,通過每個水分子的Voronoi單胞體積來區(qū)分水分子屬于液態(tài)還是屬于氣態(tài).

我們可以計算得到所有分子的Voronoi單胞體積,然后統(tǒng)計純氣態(tài)相和純液態(tài)相中分子的Voronoi單胞體積的分布,得到區(qū)分氣態(tài)和液態(tài)的一個臨界值.通過此值可以判斷相變過程中每個分子的狀態(tài)(氣態(tài)還是液態(tài)),分子的Voronoi單胞體積大于此值,認(rèn)為此分子是氣態(tài)分子,小于此值認(rèn)為是液態(tài)分子.

當(dāng)兩個分子共有同一個Voronoi面時,認(rèn)為這兩個分子是近鄰分子.當(dāng)兩個近鄰分子處于同一相(氣相或液相)時,認(rèn)為這兩個分子屬于同一個聚集體.由Voronoi方法可以給出每個分子的近鄰分子數(shù),以及兩個近鄰分子的界面面積等值.通過此方法我們可以獲得體系中液相聚集體(液滴)、氣相聚集體(氣泡)的數(shù)目,每一個液滴、氣泡的體積大小和所包含的分子數(shù)等微觀信息.

2.4 形狀因子

為了定量地描述相變過程中氣液兩相界面的形貌特征,我們計算了序參量形狀因子[33]

其中A是液相和氣相界面的總面積,A0l(A0g)是與體系中所有液滴(氣泡)總體積相等的球體的表面積.因此,為了盡可能準(zhǔn)確地估計聚集體的形狀,我們需要盡量精確地計算兩相分子所占的體積及其界面面積.雖然Voronoi單胞方法可以得到兩相分子的體積及其界面的面積,但是此界面是由多個二維平面拼接而成,形成的不是一個光滑的表面,而是一個鋸齒狀的粗糙表面,這樣求得的界面面積會偏大.

因此我們使用另一種方法來計算界面面積.液滴是由多個液態(tài)水分子堆積而成的凝聚體,計算液滴的表面積就類似于計算蛋白質(zhì)等物質(zhì)的溶液可接觸面積.把液滴中的水分子視作質(zhì)點,用兩個半徑分別為r1和r2(r2＞r1)的探測小球,在液滴表面滾動,兩個探測小球的質(zhì)心經(jīng)過的區(qū)域會形成一個封閉的厚度為d=r2-r1的殼層結(jié)構(gòu),這個殼層結(jié)構(gòu)就是兩相的界面區(qū)域.界面區(qū)覆蓋的面積依賴于兩個探測小球的半徑大小,探測小球半徑如果太小,可能會掉入液滴內(nèi)水分子的縫隙中,使得界面區(qū)比較粗糙;太大則不能很好地描述界面的形狀.通常r1應(yīng)該比水分子半徑稍大,近似為液滴內(nèi)分子的平均距離尺度.厚度d應(yīng)該遠小于液滴的外表面曲率半徑,此時界面區(qū)與液滴的接觸面積和與氣泡的接觸面積近似相等.因此只要求得界面區(qū)的體積Vinter,就可以近似得到界面區(qū)覆蓋在液滴上的面積A=Vinter/d.

我們均勻分割模擬體系來計算體系中氣相區(qū)、液相區(qū)和界面區(qū)的體積.把模擬盒子等分成許多個小網(wǎng)格,判斷每個格點處于氣態(tài)區(qū)、液態(tài)區(qū)還是界面區(qū),根據(jù)每個區(qū)域的格點數(shù)目來計算其體積.網(wǎng)格尺寸相對厚度d較小時,體積的計算可以非常精確.考慮到計算量的大小,選定劃分的網(wǎng)格的數(shù)目,根據(jù)網(wǎng)格大小給定界面區(qū)厚度d.

用每個網(wǎng)格的中心點坐標(biāo)表示此網(wǎng)格的位置,可以計算出每個格點與每個液態(tài)水分子的距離,并求出每個格點與水分子的最小距離,當(dāng)其最小距離小于r1時,此格點處于液態(tài)區(qū)域,當(dāng)其大于r2時,此格點處于氣態(tài)區(qū)域,介于兩者之間表明格點處于界面區(qū).我們給出了界面區(qū)面積以及形狀因子對r1的變化曲線,如圖2所示.

從圖2中可以看出,隨著r1的增大,界面區(qū)覆蓋在液滴上的面積先快速減小,后線性增加,而形狀因子的減小趨勢逐漸變緩.當(dāng)r1很小時,處于液態(tài)區(qū)域的點可能被判斷為處于界面區(qū)或者氣態(tài)區(qū),液滴內(nèi)存在一些空洞和凹陷,導(dǎo)致界面區(qū)面積被過度估計.隨著r1的增大,這些空洞和凹陷消失,界面區(qū)域開始變得光滑.進一步增加r1導(dǎo)致界面面積線性增加.此時的r1能較好地給出界面特征,形狀因子也處于緩慢減小的區(qū)域.

圖2 (網(wǎng)刊彩色)形狀因子(a)和界面區(qū)覆蓋液滴的面積(b)隨參數(shù)r1的變化Fig.2.(color on line)Shape factor(a)and interface area(b)as a function of param eter r1.

本文中,把模擬體系的x,y,z三個維度均等分成101份,即把整個模擬盒子均分成約1013≈106個小網(wǎng)格,每個網(wǎng)格遠小于水分子的體積.界面區(qū)的厚度與小網(wǎng)格的邊長相等,可以保證界面區(qū)被網(wǎng)格完全充滿.r1的值在界面面積線性增加的起點附近,為4.3?,比水分子尺寸稍大,符合預(yù)期.

2.5 模擬細(xì)節(jié)

本文全部模擬工作由LAMMPS分子動力學(xué)模擬軟件[39]實現(xiàn)的gNPT-MD方法完成.選用全原子水模型TIP4P/2005[40,41],體系中共包含4000個水分子.通過PPPM方法[42]來計算長程相互作用,庫侖相互作用和Lennard-Jones相互作用的截斷半徑均為9.0?.當(dāng)原子移動超過1.0?時,更新近鄰列表.通過Nosé-Hoover溫度控制器和壓強控制器[43,44]控制溫度和壓強,選用了30,55,135 atm三個壓強值.模擬的時間步長選用2 fs,x,y,z三個方向均采用周期性邊界條件.

從水的標(biāo)準(zhǔn)液態(tài)構(gòu)象出發(fā),在gNPT系綜中模擬持續(xù)升溫過程(減少β0或增加H0). 高壓(135 atm)下,每個溫度點模擬8 ns,低壓(30 atm, 55 atm)下每個溫度點模擬14 ns.用水的標(biāo)準(zhǔn)氣態(tài)構(gòu)象作為初始構(gòu)象進行降溫過程的模擬,同樣依次模擬各個溫度點,每個溫度點的模擬時長與升溫過程相同.均取最后2 ns的數(shù)據(jù)進行分析.

3 結(jié)果與討論

gNPT-MD能夠有效地探測全原子水模型的氣液相變過程,抽樣到熱力學(xué)不穩(wěn)定的共存區(qū)構(gòu)象.圖3給出了不同壓強下氣液相變的T-H曲線,每一個壓強下都有兩條模擬曲線,分別為氣態(tài)到液態(tài)的相變(降溫)和液態(tài)到氣態(tài)的相變(升溫).我們通過gNPT系綜得到的結(jié)果顯示,升溫、降溫兩條T-H曲線重合,即正則系綜下氣液相變過程中的熱滯現(xiàn)象消失,表明在gNPT系綜中,體系能連續(xù)地在兩相間轉(zhuǎn)變.曲線共分為三個區(qū)域：曲線左側(cè),溫度隨著焓值的增大而升高,此時體系處于純液相;曲線右側(cè),溫度隨著焓值的增大同樣逐漸升高,此時體系處于純氣相;曲線中間部分,溫度隨著焓值的增大而降低,此時體系處于兩相共存狀態(tài).從而,不僅是純氣體、液體相、過飽和區(qū)、spinodal點以及相共存區(qū)的微觀構(gòu)象都能在gNPT系綜中充分給出.圖中黑色虛線是30 atm下NPT系綜中的相變曲線,通過持續(xù)升溫依次模擬各溫度點得到液相到氣相的轉(zhuǎn)變(向右箭頭),通過持續(xù)降溫得到氣相到液相的轉(zhuǎn)變(向左箭頭),兩條曲線不重合,表示這里存在熱滯現(xiàn)象.在純氣態(tài)區(qū)和純液態(tài)區(qū),NPT系綜的結(jié)果與gNPT系綜的結(jié)果相符,但進入過飽和區(qū)后,非?？焖俚赝ㄟ^共存區(qū)而不能充分給出兩相共存的構(gòu)象.而且由于熱滯的存在,這種持續(xù)升降溫的NPT系綜模擬不能準(zhǔn)確地給出兩相共存溫度等相變的熱力學(xué)性質(zhì).

圖3 (網(wǎng)刊彩色)不同壓強下溫度T隨平均焓值的變化 空心圓圈表示氣態(tài)到液態(tài)的相變過程,實心星號表示液態(tài)到氣態(tài)的相變過程.彩色實線表示gNPT系綜給出的結(jié)果,帶“×”號的黑色虛線表示30 atm下NPT系綜模擬結(jié)果,向右箭頭表示其升溫過程,向左箭頭表示其降溫過程.右側(cè)給出了30 atm下相變過程中的構(gòu)象圖,紅色小球表示液態(tài)水分子,灰色小球表示氣態(tài)水分子Fig.3.(color on line)The temperatu re as a function of average enthalpy under d iff erent p ressu res:Hollow circles correspond to transition process from gas to liquid,and solid asterisks correspond to the opposite process,the colourfu l solid lines correspond to the resu lts of gNPT ensem b le,the b lack dashed lines with“×”rep resent the resu lts of NPT ensem b le at 30 atm,and the right arrow corresponds to the resu lt of heating process and the left arrow corresponds to the resu lt of cooling process.The con figurations of phase transition at 30 atm are shown on the right,the red balls rep resent liquid-like m olecu les,and the gray balls rep resent gas-likem olecules.

從圖3中可以看出,低壓時,一階相變特征明顯,氣液兩相的spinodal點的焓值相差較大,當(dāng)壓強接近TIP4P/2005水模型的臨界壓強時,一階相變特征趨于消失,兩相過渡更加光滑.隨著壓強的增大,水能達到的過熱溫度逐漸升高,這是因為在高壓時,液態(tài)水克服壓強效應(yīng)變成氣態(tài),所需要的能量更多.

圖3右側(cè)給出了30 atm時從氣態(tài)到液態(tài)的相變過程中的構(gòu)象圖.A點處于純氣相區(qū)域,此時體系中所有分子均為氣相分子,用灰色小球表示.隨著溫度降低,從純相區(qū)經(jīng)過過飽和區(qū)到達spinodal點(B點)附近時,體系中開始出現(xiàn)液相分子,用紅色小球表示.在氣液共存區(qū),液態(tài)水分子聚集成多個液滴,隨著H值的下降,液滴融合變大,浸在連通的氣泡中,如C,D兩點.在D點,體系中的氣相和液相互相分離,幾乎所有液滴融合成一個大液滴.在液態(tài)的spinodal點(E點)時,整個體系基本處于液態(tài),其中包含有多個小氣泡.F點處于純液態(tài)區(qū),所有水分子均為液相分子,然而構(gòu)象中仍有幾個氣態(tài)分子存在,可能是Voronoi單胞方法區(qū)分氣相和液相時產(chǎn)生的誤差.

在D,E兩點之間存在著從單個連通的氣泡到多個離散的小氣泡的急劇變化,這種變化隨著壓強的降低而更加顯著.這兩種構(gòu)象的焓值比較接近,但拓?fù)浣Y(jié)構(gòu)卻存在很大差別,表明低壓時,這里似乎對應(yīng)一個結(jié)構(gòu)的轉(zhuǎn)變.

為了定量地描述相變過程中液滴、氣泡數(shù)目、大小、形狀和拓?fù)浣Y(jié)構(gòu)等,基于Voronoi單胞分析方法,我們計算了如下幾個參數(shù),如圖4所示.最上方的曲線是相變區(qū)βs(H)曲線,中間的部分給出了較大的液滴、氣泡的數(shù)目(水分子數(shù)大于3的聚集體),最下方的曲線給出了最大的液滴、氣泡所含的水分子的數(shù)目.

在氣態(tài)到液態(tài)的相變過程(即焓值減小的過程)中,在氣態(tài)的過飽和區(qū),spinodal點附近,液滴的數(shù)目開始增多,而且壓強越大,液滴數(shù)目越多.進入共存區(qū)后,隨著焓值的減小,液滴的數(shù)目逐漸減少,最后全部液滴連通成一個.在液態(tài)spinodal點附近,直到液態(tài)過飽和區(qū),氣相分子不再連通,以多個小氣泡的形式出現(xiàn).圖中的紅色虛線與倒溫度曲線圍成的上下兩塊區(qū)域的面積相等,虛線對應(yīng)的溫度就是正則系綜下的共存溫度.虛線與βs(H)曲線在共存區(qū)的交點是自由能能壘的位置,近似對應(yīng)此共存溫度下的臨界核位置.

圖4 (網(wǎng)刊彩色)不同壓強下相變的βs(H)曲線(a),較大液滴(氣泡)的數(shù)目(b),最大液滴(氣泡)所含的分子數(shù)(c)隨平均焓值的變化Fig.4.(color on line)Cu rves ofβs(H)under d iff erent p ressu res(a),the num ber of bigger clusters of gas(red)and liquid(green)(b), the num ber ofm olecu les in the biggest cluster(c)as functions of average enthalpy under d iff erent p ressu res.

圖5 (網(wǎng)刊彩色)不同壓強下形狀因子在兩相共存區(qū)隨平均焓值的變化,下面的圖片給出了共存溫度點附近A,B, C三點的構(gòu)象圖Fig.5.(color on line)Shape factor dependence of average enthalpy at phase-coexisting area under different p ressures.The snapshots below show the con figu rations of A,B and C of cu rves around coexisting temperatu re.

我們計算了最大的液滴或者氣泡的形狀因子,結(jié)果如圖5所示.形狀因子值越小,越接近于1時,聚集體的形狀越接近球形;值越大,越遠離球形,聚集體表面越粗糙.A,B,C三個點分別對應(yīng)三個壓強下氣液相變過程中共存溫度處的臨界核位置,其形狀因子的值隨著壓強的升高而增大.圖下方給出了三個點的構(gòu)象圖,可以直觀地看到體系中聚集體的形狀：30 atm時,液滴非常接近球形;55 atm時,最大液滴的表面變得粗糙,開始出現(xiàn)一些分枝結(jié)構(gòu);135 atm時,最大液滴完全散落在整個體系中,全部是枝杈結(jié)構(gòu),完全遠離了球形.因此低壓時,臨界核近似于球形,近似符合CNT的假設(shè),而高壓時的狀態(tài)已經(jīng)背離了CNT的描述.

圖6 (網(wǎng)刊彩色)共存溫度(a)及表面張力隨壓強的變化(b) 圖中黃色方塊表示文獻[35]的結(jié)果,藍色菱形表示文獻[41]給出的結(jié)果,黑色星號表示文獻[45]中給出的結(jié)果,紅色圓圈表示此工作的結(jié)果.Fig.6.(color on line)Coexisting temperatu re(a)and su rface tension(b)dependence of p ressu res.The resu lts of Ref.[35](yellow square),Ref.[41](b lue d iam ond),Ref.[45](black star)are shown and the red circles correspond to the resu lts of this work.

表面張力的大小決定了臨界核對球形的偏離,我們計算了三個壓強下的共存溫度,并給出了此共存溫度處的表面張力,如圖6所示, 30,55,135 atm三個壓強下的共存溫度分別為534.0,572.0,634.0 K.文獻[35,41,45]中給出了TIP4P/2005水模型在NVT系綜下的臨界點附近的溫度變化,并根據(jù)壓強張量得到了相應(yīng)的系統(tǒng)壓強,文獻[45]中還給出了包括2740個水分子的體系的表面張力隨溫度、壓強的變化.我們模擬得到的共存溫度隨壓強的變化與文獻的結(jié)果基本相符,計算得到的表面張力也近似符合文獻中的結(jié)果.我們的結(jié)果表明高壓時氣態(tài)和液態(tài)的密度接近,氣液兩相的表面張力很小,液滴的形狀因漲落會出現(xiàn)許多枝杈結(jié)構(gòu).低壓時氣態(tài)和液態(tài)的密度相差較大,表面張力較大,液滴近似為球形.

4 結(jié) 論

為了研究氣液相變過程中氣液共存狀態(tài)下的微觀細(xì)節(jié),我們通過把倒溫度改為焓依賴的函數(shù)的方式,給出了廣義的等溫-等壓系綜,然后在分子動力學(xué)標(biāo)準(zhǔn)軟件包中實現(xiàn)了gNPT-MD方法.此方法不需要使用副本交換等技術(shù)就可以模擬復(fù)雜的大系統(tǒng)的相變過程,能有效地抽樣得到熱力學(xué)不穩(wěn)定的共存區(qū)的構(gòu)象.基于此方法,我們模擬了全原子水模型TIP4P/2005在不同壓強下的氣液相變過程.

可以發(fā)現(xiàn),用gNPT-MD得到的純氣態(tài)和純液態(tài)區(qū)的相變曲線能與常規(guī)NPT系綜下得到結(jié)果很好地符合,而且gNPT-MD方法的結(jié)果消除了NPT系綜中的過飽和熱滯現(xiàn)象,能夠給出穩(wěn)定存在的共存區(qū)的構(gòu)象.低壓時,相變曲線具有明顯的一階相變的特征,隨著壓強的升高,在接近模型的臨界壓強的高壓時,一階相變特征逐漸消失,相變開始具有二階相變特征.

通過Voronoi單胞分析方法,我們分析了相變過程中形成的氣泡、液滴的變化,發(fā)現(xiàn)在氣態(tài)spinodal點附近,液滴的數(shù)目隨著壓強的升高而增多.我們改進了計算氣液兩相界面面積的方法,通過計算液滴的表面可達面積的方法,得到了更加精確的面積值.計算了序參量形狀因子fs,來描述體系中聚集體的形貌特征,發(fā)現(xiàn)共存區(qū)的微觀結(jié)構(gòu)在不同壓強下存在著明顯的差別.本文的結(jié)果表明：低壓時氣液兩相的密度相差較大,因此表面張力較大,臨界核能被維持成一個有限的球形,此時的形狀因子值較小,更接近1;高壓時氣態(tài)和液態(tài)的spinodal點非常接近,氣液兩相界面的表面張力很小,臨界核會出現(xiàn)許多枝杈結(jié)構(gòu),遠離了球形,此時的形狀因子值較大.

感謝中國科學(xué)院物理研究所楊成博士及中國科學(xué)院大學(xué)物理學(xué)院張傳彪的討論與幫助.

[1]Erdem ir D,Lee A Y 2009 Acc.Chem.Res.42 621

[2]Sleu telM,Lu tsko J,van D riessche A E,Du rán-O livencia M A,M aes D 2014 Nat.Comm un.5 5598

[3]Auer S,Frenkel D 2004 Annu.Rev.Phys.Chem.55 333

[4]Toxvaerd S 2015 J.Chem.Phys.143 154705

[5]Debenedetti P G 2006 Nature 441 168

[6]Gasser U,W eeks E R,Schofield A,Pusey P N,Weitz D A 2001 Science 292 258

[7]Yasuoka K,M atsum oto M 1998 J.Chem.Phys.109 8451

[8]Yasuoka K,M atsum oto M 1998 J.Chem.Phys.109 8463

[9]M yerson A S,Trout B L 2013 Science 341 855

[10]Savage J R,D insm ore A D 2009 Phys.Rev.Lett.102 198302

[11]Sleu telM,van D riessche A E 2014 Proc.Natl.Acad.Sci. 111 E546

[12]de Yoreo J 2013 Nature M ater.12 284

[13]Yarom M,M arm ur A 2015 Adv.Colloid In terface Sci. 222 743

[14]Du?ka M,Něm ec T,H rubyJ,V in?V,P lankováB 2015 EPJW eb Conf.92 02013

[15]Schenter G K,K athm ann SM,Garrett B C 1999 Phys. Rev.Lett.82 3484

[16]Reguera D,Reiss H 2004 Phys.Rev.Lett.93 165701

[17]Bhim alapuram P,Chakrabarty S,Bagchi B 2007 Phys. Rev.Lett.98 206104

[18]Rane K S,M u rali S,Errington J R 2013 J.Chem.Theory Com put.9 2552

[19]PlankováB,Vin?V,H rubyJ,Du?ka M,Něm ec T,CelnyD 2015 EPJW eb Conf.92 02071

[20]M cG rath M J,K uo I F W,Ghogom u J N,M undy C J, Siepm ann J I 2011 J.Phys.Chem.B 105 11688

[21]M alolepsza E,K im J,K eyes T 2015 Phys.Rev.Lett. 114 170601

[22]Kuo IF W,M undy C J 2004 Science 303 658

[23]Nagata Y,Usui K,Bonn M 2015 Phys.Rev.Lett.115 236102

[24]Zahn D 2004 Phys.Rev.Lett.93 227801

[25]Panagiotopou los A Z 1987 M ol.Phys.61 813

[26]M ou?ka F,Nezbeda I 2013 F luid Phase Equilib.360 472

[27]Trejos V M,Gil-V illegas A,M artinez A 2013 J.Chem. Phys.139 184505

[28]Cho W J,K im J,Lee J,Keyes T,Straub J E,K im K S 2014 Phys.Rev.Lett.112 157802

[29]K im J,Keyes T,Straub J E 2010 J.Chem.Phys.132 224107

[30]M a?olepsza E,Secor M,K eyes T 2015 J.Phys.Chem. B 119 13379

[31]Lu Q,K im J,Straub J E 2013 J.Chem.Phys.138 104119

[32]Xu S,Zhou X,Ouyang Z C 2012 Comm un.Com put. Phys.12 1293

[33]Jeong S,Jho Y,Zhou X 2015 Sci.Rep.5 15955

[34]G loor G J,Jackson G,B las F J,de M iguel E 2005 J. Chem.Phys.123 134703

[35]Vega C,de M iguel E 2007 J.Chem.Phys.126 154707

[36]K um ar V S,K um aran V 2005 J.Chem.Phys.123 114501

[37]Zhu H X,Thorpe S M,W ind le A H 2001 Philos.M ag. A 81 2765

[38]Oger L,Gervois A,Troadec J P,Rivier N 1996 Philos. M ag.B 74 177

[39]P lim p ton S 1995 J.Com put.Phys.117 1

[40]Abascal J L,Vega C 2005 J.Chem.Phys.123 234505

[41]Vega C,Abascal J L F,Nezbeda I 2006 J.Chem.Phys. 125 034503

[42]Beckers J V L,Lowe C P,de Leeuw S W 1998 M ol. Sim u l.20 369

[43]NoséS 1984 J.Chem.Phys.81 511

[44]Hoover W G 1985 Phys.Rev.A 31 1695

[45]Alejand re J,Chapela G A 2010 J.Chem.Phys.132 014701

(Received 29 March 2017;revised manuscript received 4 May 2017)

Vapor-liquid coexisting morphology of all-atom water model through generalized isothermal isobaric ensemble molecular dynamics simulation?

Yin Ling-Kang1)Xu Shun2)Seongm in Jeong3)Yongseok Jho3)Wang Jian-Jun4)Zhou Xin1)?

1)(School of Physical Sciences,University of Chinese Academ y of Sciences,Beijing 100049,China)
2)(Com pu ter Network Inform ation Center,Chinese Academ y of Sciences,Beijing 100190,China)
3)(Center for Soft and Living M atter,Institute for Basic Science,IBS,U lsan 44919,Repub lic of Korea)
4)(Institu te of Chem istry,Chinese Academ y of Sciences,Beijing 100190,China)

Exp loring the atom-scale details such asmorphology of coexisting phase during phase transitions is very im portant for understanding their microscopic m echanism.While most theories,such as the classic nucleation theory,usually over-sim p lify the character of the critical nucleus,like the shape,structure,and most current experiment techniques are hard ly to cap ture the instantaneousmicroscopic details,the atomistic molecular dynamics(MD)or Monte Carlo(MC) simulation provides a promise to detect the intermediate process of phase transitions.However,the standard canonicalensemble MD/MC simulation technique can not sufficiently sample the instantaneous(unstable in therm odynam ics) coexistent phase.Therefore,the MC in the general canonical ensemb le,such as general isothermal-volume ensemble (gNVT),combined with the enhanced sampling techniques,such as the rep lica exchange(RE)method,was presented to stabilize then to su ffi ciently sam p le the atom ic con formations of the phase coexistence.Due to the lim it of the RE, the RE-MC simulation on gNVT is usually app lied in sm aller system s.In this paper,we fi rst extend the gNVT-based MC simulation to the MD in the generalized isothermal-isobaric ensemble(gNPT)and very simplyimplement it in the standard atomic MD soft packages withoutm odifying the code,so that we can use these packages in MD simulation of realistic systems.Then we simulate the vapour-liquid phase transition of all-atom ic water model.At least at not very low pressures,we find that the individual gNPT simulation is already enough to reach equilibrium in any region of the phase transition,not only in the normal liquid and vapour regions,but in the super-saturation regions,and even in the vapour-liquid coexistent regions.The obtained energy-temperature curve in the cooling gNPT wellm atches with that in the heating procedure without any hysteresis.It indicates that it is not necessary to use the RE technique in the gNPT,and the interm ediate states during phase transitions in larger system s can be effectively simulated by a series of independent individual gNPT-MD simulations in the standard soft packages.We also p ropose a method to accurately determ ine the interface between the two phases in the coexistence,then provide a quantitativem easurement about the interface tension and themorphology of the coexistent phase in the larger all-atomic water at various temperatures and pressures.The results show that the liquid droplet(or vapour bubble)at the low pressure is close to a spheredue to the larger interface tension,as expectation of the classic nucleation theory of the first-order phase phase transition,but becom esm ore and m ore irregu lar as the decrease of the interfacial tension as increasing the pressure to approach to the critical p ressure,where the phase transition is the second order one.

generalized canonical ensemble,gas liquid transition,gas liquid coexistence,molecular dynam ics

PACS:61.20.Ja,64.60.Q-,64.70.F-,31.15.xv DO I:10.7498/aps.66.136102

?國家自然科學(xué)基金(批準(zhǔn)號:11574310)資助的課題.

?通信作者.E-m ail:xzhou@ucas.ac.cn

PACS:61.20.Ja,64.60.Q-,64.70.F-,31.15.xv DO I:10.7498/aps.66.136102

*Pro ject supported by the National Natural Science Foundation of China(G rant No.11574310).

?Corresponding author.E-m ail:xzhou@ucas.ac.cn