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

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

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

2)(中國科學院計算機網絡信息中心,北京 100190)

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

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

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

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

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

2)(中國科學院計算機網絡信息中心,北京 100190)

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

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

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

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

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

1 引 言

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

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

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

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

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

2 研究方法

2.1 GCE與gNPT系綜

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

NVT系綜的構象分布函數可以寫為

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

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

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

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

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

這里P和V分別是外部壓強和系統體積,H0是常數.同樣可以有gNPT系綜中焓依賴的熱源倒溫度：

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

2.2 gNPT系綜與NPT系綜的關系

我們以gNPT與NPT系綜為例說明廣義系綜與正則系綜的關系,GCE和NVT系綜的關系基本類似.任意系統的構象空間性質由其態密度eS(H),或者內稟倒溫度函數

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

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

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

圖1 (網刊彩色)典型的相變過程的倒溫度曲線β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).

由此我們容易得到任意溫度T下NPT系綜的H空間分布,PNPT(H)∝e-H/T+S(H),以及自由能,

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

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

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

2.3 Voronoi單胞分析方法

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

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

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

2.4 形狀因子

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

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

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

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

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

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

圖2 (網刊彩色)形狀因子(a)和界面區覆蓋液滴的面積(b)隨參數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個小網格,每個網格遠小于水分子的體積.界面區的厚度與小網格的邊長相等,可以保證界面區被網格完全充滿.r1的值在界面面積線性增加的起點附近,為4.3?,比水分子尺寸稍大,符合預期.

2.5 模擬細節

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

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

3 結果與討論

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

圖3 (網刊彩色)不同壓強下溫度T隨平均焓值的變化 空心圓圈表示氣態到液態的相變過程,實心星號表示液態到氣態的相變過程.彩色實線表示gNPT系綜給出的結果,帶“×”號的黑色虛線表示30 atm下NPT系綜模擬結果,向右箭頭表示其升溫過程,向左箭頭表示其降溫過程.右側給出了30 atm下相變過程中的構象圖,紅色小球表示液態水分子,灰色小球表示氣態水分子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點的焓值相差較大,當壓強接近TIP4P/2005水模型的臨界壓強時,一階相變特征趨于消失,兩相過渡更加光滑.隨著壓強的增大,水能達到的過熱溫度逐漸升高,這是因為在高壓時,液態水克服壓強效應變成氣態,所需要的能量更多.

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

在D,E兩點之間存在著從單個連通的氣泡到多個離散的小氣泡的急劇變化,這種變化隨著壓強的降低而更加顯著.這兩種構象的焓值比較接近,但拓撲結構卻存在很大差別,表明低壓時,這里似乎對應一個結構的轉變.

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

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

圖4 (網刊彩色)不同壓強下相變的βs(H)曲線(a),較大液滴(氣泡)的數目(b),最大液滴(氣泡)所含的分子數(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 (網刊彩色)不同壓強下形狀因子在兩相共存區隨平均焓值的變化,下面的圖片給出了共存溫度點附近A,B, C三點的構象圖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.

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

圖6 (網刊彩色)共存溫度(a)及表面張力隨壓強的變化(b) 圖中黃色方塊表示文獻[35]的結果,藍色菱形表示文獻[41]給出的結果,黑色星號表示文獻[45]中給出的結果,紅色圓圈表示此工作的結果.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系綜下的臨界點附近的溫度變化,并根據壓強張量得到了相應的系統壓強,文獻[45]中還給出了包括2740個水分子的體系的表面張力隨溫度、壓強的變化.我們模擬得到的共存溫度隨壓強的變化與文獻的結果基本相符,計算得到的表面張力也近似符合文獻中的結果.我們的結果表明高壓時氣態和液態的密度接近,氣液兩相的表面張力很小,液滴的形狀因漲落會出現許多枝杈結構.低壓時氣態和液態的密度相差較大,表面張力較大,液滴近似為球形.

4 結 論

為了研究氣液相變過程中氣液共存狀態下的微觀細節,我們通過把倒溫度改為焓依賴的函數的方式,給出了廣義的等溫-等壓系綜,然后在分子動力學標準軟件包中實現了gNPT-MD方法.此方法不需要使用副本交換等技術就可以模擬復雜的大系統的相變過程,能有效地抽樣得到熱力學不穩定的共存區的構象.基于此方法,我們模擬了全原子水模型TIP4P/2005在不同壓強下的氣液相變過程.

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

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

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

[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

?國家自然科學基金(批準號: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