分子模擬教程徑向分布函數(shù)

1、第三章 分子模擬方法,1. 蒙特卡羅(Monte Carlo)方法基礎(chǔ) 2. 分子動力學(Molecular Dynamics)方法基礎(chǔ) 3. 程序講解,本章具體講解內(nèi)容,掌握分子模擬方法的必備知識：,編程技能 (Fortran or C/C++) 統(tǒng)計物理學（統(tǒng)計力學）： 統(tǒng)計物理學基礎(chǔ)； 系綜原理； 非平衡統(tǒng)計力學基礎(chǔ)； 漲落理論分子熱力學 ：分子間相互作用理論； 分布函數(shù)理論氣體分子運動論其它,分

2、子模擬的目的:,為什么要進行分子模擬？,將分子聚集體的性質(zhì)與如下方面相聯(lián)系：分子的微觀相互作用分子聚集體的結(jié)構(gòu)分子的動力學過程,分子模擬對實驗進行補充，使我們能夠：預測現(xiàn)有或新材料的性質(zhì)在分子水平研究宏觀現(xiàn)象獲得實驗無法或難以發(fā)現(xiàn)的東西,什么是計算機分子模擬方法？,分子模擬的定義： 統(tǒng)計力學基本原理出發(fā)，將一定數(shù)量的分子輸入計算機內(nèi)進行分子微觀結(jié)構(gòu)的測定和宏觀性質(zhì)的計算。,按照獲得微觀態(tài)的方法不同，分子模擬分為：蒙特

3、卡羅方法 (Monte Carlo, MC)分子動力學方法 (Molecular Dynamics, MD) (3) 混合方法 (hybrid method，HM),計算機分子模擬的發(fā)展歷史：,1. 蒙特卡羅方法（MC）1953 Metropolis, Ulam, Rosenbluth and Tell， Los Alamos National LabMont

4、e Carlo simulation of hard sphere.2. 分子動力學方法（MD）1957Alder and Wainwrigth，Livermore LabMolecular dynamics simulation of hard spheres.,微觀與宏觀,分子模擬在微觀尺度與實驗室的宏觀世界之間起著橋梁的作用：,給定分子間的相互作用 “準確”預測研究體系的性質(zhì),MC與MD的區(qū)別：,MC:

5、 構(gòu)型平均，不包含動力學部分； 利用概率行走產(chǎn)生微觀態(tài)。MD:時間平均，產(chǎn)生動力學性質(zhì)；利用運動軌線隨時間的變化來產(chǎn)生一系列微觀態(tài)。,計算機分子模擬的發(fā)展歷史（續(xù)）：,從上個世紀九十年代初期以來，計算機模擬技術(shù)得到了飛速發(fā)展，主要基于三個方面的發(fā)展：分子力場的發(fā)展（基石） （Amber，OPLS、Compass）原子間的鍵長、鍵角、分子間的內(nèi)聚能等模擬算法（途徑）計算機硬件（工具）,HPCx,計算

6、機分子模擬的特點：,原子水平的模擬計算機實驗檢驗理論、篩選實驗科學研究中的第三種方法,分子模擬中涉及的幾個基本概念：,模擬計算盒子或模擬胞腔,Simulation box (cell),裝有一定數(shù)目流體分子的研究對象，它是我們要研究的宏觀體系的縮微模型。,立方形胞腔,,,,周期性邊界條件(Periodic boundary condition, PBC),在小體系中，邊界效應(yīng)總是很顯著。在包含1000個原子的簡單立方晶體中

7、－488個原子處于邊界上。在包含1000000個原子的簡單立方晶體中－仍然有 6%的原子在邊界上。,在模擬中，考慮具有真實邊界的對象，不切合實際： 增強了有限尺寸效應(yīng) 人為造成的邊界會影響流體的性質(zhì),當某個粒子運動出模擬盒子的某一邊界時，另外一個影像粒子從另一對立邊界進入到此盒子中。,周期性邊界條件(Periodic boundary condition, PBC),本體體系的近似：中心盒子在X，Y和Z方向無限擴展；消除人為形成

8、邊界的表面效應(yīng)；保證中心盒子中的粒子數(shù)恒定。只需要跟蹤中心盒子中各粒子的運動。,周期性邊界條件的算法：,采用數(shù)學函數(shù)：,FLOOR(r/L): 返回不超過r/L的最大整數(shù),FLOOR (4.8) has the value 4.FLOOR (-5.6) has the value -6.,采用數(shù)學函數(shù)：,r/L>0, ANINT(r/L) = AINT(r/L+0.5),r/L?0, ANINT(r/L) = AINT(

9、r/L-0.5),周期性邊界條件的算法：,y,最小影像轉(zhuǎn)化原理(Minimum image convention),定義：,中心元胞中的一個粒子只與此元胞中的其它N－1個粒子，或它們的最近鄰影像發(fā)生相互作用。,適用條件：,粒子間相互作用勢能的截斷距離必須不大于模擬中心元胞長度的一半。,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,此兩粒子與中心粒子的距離相等，但是：黑色球發(fā)生作用綠色球不發(fā)生作用,,此兩粒子是與中心

10、原子相互作用的最近鄰影像,最小影像轉(zhuǎn)化原理的算法：,采用數(shù)學函數(shù)：,r/L>0, ANINT(r/L) = AINT(r/L+0.5),r/L?0, ANINT(r/L) = AINT(r/L-0.5),截斷勢能(Truncating the Potential),,,本體體系采用周期性邊界條件描述：不可能將所有粒子與它們影像粒子間的相互作用全都計算。 必須在不大于中心盒子長度的一半處進行截斷，以便與最小影像轉(zhuǎn)化原理一致。

11、粒子間的相互作用主要來自于截斷范圍內(nèi)，而范圍外的貢獻很小，可忽略不計。,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,截斷范圍內(nèi)的相互作用,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,截斷勢能函數(shù)的形式：,簡單截斷勢能函數(shù)(Truncated Potential)：,缺點：,rc: 截斷距離或半徑,勢能在截斷處不連續(xù),當一對分子穿越邊界時，總能量不守恒。 分子間力在截斷處為無窮大，MD運

12、動過程不穩(wěn)定。,,忽略截斷半徑之外的所有作用,位移截斷勢能函數(shù)(Shifted and Truncated Potential)：,,缺點：,分子間力仍然在截斷處不連續(xù)。,優(yōu)點：,勢能在截斷處連續(xù)，但不影響分子間力的大小 分子間力在截斷處不為無窮大,截斷勢能函數(shù)的形式：,常用于MC和MD模擬中,,,,位移－力截斷勢能函數(shù)(Shifted-Force Potential)：,常用于MD模擬中,優(yōu)點：,勢能和分子間力均在截斷處連續(xù),截斷勢

13、能函數(shù)的形式：,,截斷勢能函數(shù)的對比：,位移－力截斷勢能,簡單截斷勢能函數(shù),一、Monte Carlo模擬方法基礎(chǔ)：,亦稱統(tǒng)計模擬或隨機抽樣方法，statistical simulation method ?利用隨機數(shù)進行數(shù)值模擬的方法,Monte Carlo名字的由來：,是由Metropolis在二次世界大戰(zhàn)期間提出的：Manhattan計劃，研究與原子彈有關(guān)的中子輸運過程；,Nicholas Metropolis (1915-19

14、99),Monte-Carlo, Monaco,投硬幣，擲骰子,Monte Carlo方法計算Pi值,隨機數(shù)的定義和特性,什么是隨機數(shù)？,單個的數(shù)字不是隨機數(shù);,是指一個數(shù)列，其中的每一個體稱為隨機數(shù)，其值與數(shù)列中的其它數(shù)無關(guān)；,在一個均勻分布的隨機數(shù)中，每一個體出現(xiàn)的概率是均等的；,例如：在[0,1]區(qū)間上均勻分布的隨機數(shù)序列中，0.00001與0.5出現(xiàn)的機會均等,隨機數(shù)應(yīng)具有的基本特性,隨機數(shù)序列應(yīng)是獨立的、互不相關(guān)的(uncor

15、related)：,即序列中的任一子序列應(yīng)與其它的子序列無關(guān)；,長的周期(long period)：,均勻分布的隨機數(shù)應(yīng)滿足均勻性(Uniformity)：,隨機數(shù)序列應(yīng)是均勻的、無偏的，即：如果兩個子區(qū)間的“面積”相等，則落于這兩個子區(qū)間內(nèi)的隨機數(shù)的個數(shù)應(yīng)相等。,例如：對[0,1)區(qū)間均勻分布的隨機數(shù)，如果產(chǎn)生了足夠多的隨機數(shù)，而有一半的隨機數(shù)落于區(qū)間[0,0.1]?不滿足均勻性,如果均勻性不滿足，則會出現(xiàn)序列中的多組隨機數(shù)相關(guān)的情況

16、?均勻性與互不相關(guān)的特性是有聯(lián)系的,實際應(yīng)用中，隨機數(shù)都是用數(shù)學方法計算出來的，這些算法具有周期性，即當序列達到一定長度后會重復；,有效性（Efficiency):,模擬結(jié)果可靠,?模擬產(chǎn)生的樣本容量大,?所需的隨機數(shù)的數(shù)量大,?隨機數(shù)的產(chǎn)生必須快速、有效，最好能夠進行并行計算。,隨機數(shù)與隨機數(shù)發(fā)生器,得到一個可能的隨機數(shù)序列，是在計算機上實現(xiàn)Monte Carlo方法的關(guān)鍵,隨機數(shù)的產(chǎn)生方法：,[0,1]區(qū)間上均勻分布的隨機數(shù)是Mon

17、te Carlo模擬的基礎(chǔ)，服從任意分布的隨機數(shù)序列可以用[0,1]區(qū)間均勻分布的隨機數(shù)序列作適當?shù)淖儞Q或舍選后求得。,??(0,1) ? 2 ?-1?(-1,1),利用隨機數(shù)表，如Tippett于1972年發(fā)表的隨機數(shù)表； 占用太多的計算機內(nèi)存采用物理方法，如利用電子線路的熱噪聲等； 昂貴而且不便重復,:,?偽隨機數(shù)（Pseudo-Random Number),遞推到一定次數(shù)后，出現(xiàn)周期性的重復現(xiàn)象。,利用數(shù)學遞推

18、公式,一旦公式和初值定下來，整個隨機數(shù)序列便被確定下來，而且每一個隨機數(shù)只被它前面的那個數(shù)唯一確定，因此這類隨機數(shù)并不是真正的隨機數(shù)。,Monte Carlo方法基本思想,當所求的問題是某種事件出現(xiàn)的概率，或是某個隨機變量的期望值時，它們可以通過某種“隨機試驗”的方法，得到這種事件出現(xiàn)的頻率和概率，或者得到這個隨機變量的統(tǒng)計平均值，并用它們作為問題的解。,Monte Carlo方法解決的問題,問題本身是確定性問題，要求我們?nèi)ふ乙粋€隨機

19、過程，使該隨機過程的統(tǒng)計平均就是所求問題的解。,問題本身就是隨機過程，我們可以根據(jù)問題本身的實際物理過程來進行計算機模擬和跟蹤，并采用統(tǒng)計方法求得問題的解。,Monte Carlo方法的特點,計算的收斂性和收斂速度均與問題的維數(shù)無關(guān)，適合解決高維問題。,對問題的適應(yīng)能力強。,收斂速度僅為樣本數(shù)的-1/2次，因而計算耗時大。,Monte Carlo方法的應(yīng)用舉例：,計算積分：,常用的積分方法求解：,將積分區(qū)域[a,b]均勻地劃分成N各分區(qū)

20、間，則積分結(jié)果可近似地表示成：,Δx = (b-a)/N,簡單的Monte Carlo積分方法求解：,利用均勻分布的隨機數(shù)發(fā)生器，從[a,b]區(qū)間產(chǎn)生一系列隨機數(shù)xi，i=1, 2, ..., N,其中 X為均勻分布，并且 X?[a,b],近似求解E[g(X)]:,近似求解積分:,隨機抽樣,當我們用簡單Monte Carlo計算積分時，若該函數(shù)為常數(shù)函數(shù)，g(x)=constant，則取樣數(shù)不管多少，準確度為100％。如果在積分區(qū)間

21、內(nèi)，g(x)為一平滑函數(shù)，則簡單Monte Carlo方法較為準確，反之，如果g(x)的變動很劇烈，則簡單Monte Carlo方法的誤差會變大。,說明：,,,重要性Monte Carlo抽樣方法,在 g(x) 變化劇烈時，如果以Monte Carlo方法取樣，最好依據(jù)g(x)的大小來決定取樣率。當|g(x)|的值較大時，對∫g(x)dx的貢獻也較大，如果沒被選中，則結(jié)果的誤差極大。解決方式：改變 x被選中的機率，讓|g(x)| 值

22、較大的點被選中的機率增加。采用權(quán)重分布函數(shù)(Weight distribution function) w(x) :決定每個x被選中的機率。,重要性抽樣的定義：根據(jù)一定的分布形式進行的隨機抽樣。,w(x)必須歸一化，即在積分區(qū)間內(nèi)∫w(x)dx=1。由于 x 的選取已被 w(x) 扭曲，所以計算積分時要把這部分[還]回去：若一共取樣了N個x，則積分值為：,,,重要性Monte Carlo抽樣方法,Metropolis Monte C

23、arlo方法,我們所模擬的系統(tǒng)最終要達到的平衡分布是Boltzman分布：,Boltzmann 概率分布函數(shù)：,,我們?nèi)绻軌虍a(chǎn)生這種分布，我們就能夠計算系統(tǒng)的大多數(shù)性質(zhì)，但這是不可能的，因為我們不知道Z的值，但是對于任意兩個狀態(tài)，我們有：,,可以在相空間中構(gòu)造一個馬爾科夫鏈，使相空間中的樣本點隨著鏈的增長逐步趨近于Boltzman分布。,一個序列x0, x1, x2, …,xn,如果對任何n都有：,,則此序列是一個Markov鏈。,要

24、求：,任何一次實驗的結(jié)果依賴于前一次的試驗，并且近依賴于前一次的試驗。,馬爾科夫(Markov)鏈,可以證明：通過構(gòu)造Markov鏈，體系中最終的平衡分布就是Boltzman分布,Metropolis平衡條件 (Detailed balance condition)：,平衡條件：,系統(tǒng)處于狀態(tài)X的概率正比于其Boltzman因子：,,如果?是對稱的:,Metropolis Monte Carlo方法的算法：,,,,給出一個初始狀態(tài)，并計

25、算系統(tǒng)的能量Eold,隨機產(chǎn)生一個新狀態(tài)，并計算新系統(tǒng)的能量Enew,如果ΔE＝(Enew?Eold)<0, 則接受新狀態(tài)并回到 b),如果ΔE＝(Enew?Eold)>0, 則計算Boltzman因子：,在(0, 1)區(qū)間上產(chǎn)生一個均勻分布的隨機數(shù)? ；,如果,則接受新狀態(tài)并回到b),否則保留原值并回到b),1. 正則系綜蒙特卡羅模擬方法(Canonical MC Simulation),具有確定的粒子數(shù)N、溫度T和體積V

26、,,,,對于含N個粒子的系統(tǒng)，位型（構(gòu)型）的配分函數(shù)：,某個特定構(gòu)型的發(fā)生概率為 PNVT(rN),1. 正則系綜蒙特卡羅模擬方法(Canonical MC Simulation),Monte Carlo模擬中任一物理量的計算：,位型積分,概率密度,系統(tǒng)處于位型{rN}的概率密度,,給出一個初始狀態(tài)，并計算系統(tǒng)的能量Uold,隨機產(chǎn)生一個新狀態(tài)，并計算新系統(tǒng)的能量Unew,如果ΔU＝(Unew?Uold)<0, 則接受新狀態(tài)并回到

27、 b),如果ΔU＝(Unew?Uold)>0, 則計算Boltzman因子：,在(0, 1)區(qū)間上產(chǎn)生一個均勻分布的隨機數(shù)? ；,如果,則接受新狀態(tài)并回到b),否則保留原值并回到b),正則系綜MC模擬算法的組織：,,正則系綜MC模擬算法的流程圖：,給定每個分子的初始位置，ri(0),隨機選取一個分子，并隨機移動到新的位置,計算移動前后的系統(tǒng)能量變化ΔU,拒絕移動,,,,,ΔU?0 ?,,Exp(ΔU)?? (0,1)?,,,統(tǒng)計系

28、統(tǒng)的熱力學性質(zhì)及其它物理量,,,,統(tǒng)計性質(zhì)不變？,,打印結(jié)果，結(jié)束,Yes,No,,,,No,No,,接受移動,,,,,Yes,Yes,,大約循環(huán)107到108次,Monte Carlo模擬中幾個熱力學量的計算：,N個粒子系統(tǒng)中的總勢能：,假設(shè)采用截斷勢能函數(shù)：,Uc：截斷范圍內(nèi)的總勢能；Ulrc：截斷半徑外對勢能的長程校正(Long-range correction),對于LJ流體：,,,含N個粒子系統(tǒng)中的壓力：,Wc：截斷范圍內(nèi)的

29、總維里項(Virial)；,Plrc：截斷半徑外對壓力的長程校正,,Read simulationparameters,Start,Initialize positionsof all particles,New simulation?,Read oldconfiguration,Monte Carloloop,Stop,yes,no,Monte Carlo loop Subroutine,Start,Stop,Trial m

30、ove,Satisfy Metropolisrule?,Accept thetrial move,Update energyand virial,Sample thepressure,End ofsimulation?,yes,no,,yes,no,Main program,正則系綜MC模擬程序基本結(jié)構(gòu)：,正則系綜MC模擬程序F11講解（LJ, NVT）：,** READ INPUT DATA **,初始狀態(tài)：,READ (

31、*,‘(A)’) Title ! 運行作業(yè)題目READ (*,*) NStep ! 運行步數(shù)READ (*,*) Iprint ! 打印步數(shù)READ (*,*) Isave ! 保存步數(shù)READ (*,*) Iratio ! 調(diào)整步數(shù)READ (*,‘(A)’) CNFile ! 位型文件READ (*,*) Dens

32、 ! 對比密度READ (*,*) Temp ! 對比溫度 READ (*,*) Rcut ! 對比截斷半徑,,無因次量：,,正則系綜MC模擬程序F11講解（LJ, NVT）：,量綱變換：,Beta = 1.0 / TempSigma = ( Dens / Real ( N ) ) ** ( 1.0 / 3.0 )Rmin = 0.70 * Sigma

33、 ！判斷粒子發(fā)生重疊時的距離Rcut = Rcut * Sigma ！截斷半徑DRmax = 0.15 * Sigma ！隨機移動的最大距離DensLJ = DensDens = Dens / ( Sigma ** 3 )IF ( Rcut .GT. 0.5 ) Stop ' Cut-Off Too Large ',,模擬盒子的邊長為1,,,

34、** Read Initial Configuration *Call ReadCN（CNFile）,正則系綜MC模擬程序F11講解（LJ, NVT）：,初始位型：,參閱程序F23,Call FCC,需要自己給定所有粒子初始位置,面心立方 (face-centered cubic, FCC)：,正則系綜MC模擬程序講解（LJ, NVT）：,長程校正：,Sr3 = ( Sigma / Rcut ) ** 3Sr9 = Sr3 ** 3

35、Vlrc12 = 8.0 * Pi * DensLJ * Real ( N ) * Sr9 / 9.0Vlrc6 = - 8.0 * Pi * DensLJ * Real ( N ) * Sr3 / 3.0Vlrc = Vlrc12 + Vlrc6Wlrc12 = 4.0 * Vlrc12Wlrc6 = 2.0 * Vlrc6Wlrc = Wlrc12 + Wlrc6,算法：能量求和：,C

36、all Sumup ( Rcut, Rmin, Sigma, Ovrlap, V, W )If ( Ovrlap ) Stop ' Overlap In Initial Configuration 'Vs = ( V + Vlrc ) / Real ( N )Ws = ( W + Wlrc ) / Real ( N )Ps = Dens * Temp + W + WlrcPs = Ps * Sigma **

37、 3,** Check For Acceptance **DeltV = Vnew - VoldDeltW = Wnew - WoldDeltVb = Beta * DeltVIf ( DeltVb .Lt. 75.0 ) Then If ( DeltV .Le. 0.0 ) Then V = V + DeltV W = W + DeltW

38、 RX(i) = RXInew ………….. Acatma = Acatma + 1.0 Elseif ( Exp ( - DeltVb ) .Gt. Ranf ( Dummy ) ) Then V = V + DeltV W = W + DeltW RX(i) = RXInew …………… Acatma = Acat

39、ma + 1.0 EndifEndifAcm = Acm + 1.0,算法：Metropolis算法,算法：勢能的計算,DO 100 J = 1, NIF ( I .NE. J ) ThenRXIJ = RXI - RX(J) ………RXIJ = RXIJ - Anint ( RXIJ ) ……………RIJSQ = RXIJ * RXIJ + RYIJ * RYIJ + RZIJ * RZIJIF ( R

40、IJSQ .LT. RcutSQ ) THENSR2 = SIGSQ / RIJSQSR6 = SR2 * SR2 * SR2VIJ = SR6 * ( SR6 - 1.0 )WIJ = SR6 * ( SR6 - 0.5 )V = V + VIJW = W + WIJ ENDIFENDIF100 ContinueV = 4.0 * VW = 48.0 * W / 3.0,最小影像原理,算法：Metrop

41、olis算法,RXIOld = RX(I) …………..** Calculate The Energy Of I In The Old Configuration **Call Energy ( RXIold, RYIold, RZIold, I, Rcut, Sigma, Vold, Wold )C ** Move I And Pickup The Central Image **RXInew = RXIold + ( 2

42、.0 * Ranf ( DUMMY ) - 1.0 ) * DRmax……………….. RXInew = RXInew - Anint ( RXInew )………………** Calculate The Energy Of I In The New Configuration **Call Energy ( RXInew, RYInew, RZInew, I, Rcut, Sigma,Vnew, Wne

43、w ),隨機移動,** Adjust Maximum Displacement **Ratio = Acatma / Real ( N * Iratio )If ( Ratio .Gt. 0.5 ) Then DRmax = DRmax * 1.05Else DRmax = DRmax * 0.95EndifAcatma = 0.0EndifIf ( Mod ( Step, Iprint ) .

44、Eq. 0 ) Then,算法：步長的調(diào)整,分子模擬中徑向分布函數(shù)g(r)的計算：,計算表達式：,物理意義：流體中距一個分子為 r 處出現(xiàn)另一個分子的幾率密度，它反映流體中短程有序的特點。,,,,,r,,r+?r,,,在模擬中通常在盒子長度一半的范圍內(nèi)，考察g(r)隨距離的變化。,分子模擬中徑向分布函數(shù)g(r)的算法：,第一步：計算球殼層間的距離，并初始化一些變量：,Ngr: g(r)的統(tǒng)計次數(shù),Delg: 球殼層間的距離,nhi

45、st: 球殼層的層數(shù),,模擬開始時需要給定： nhist 以及統(tǒng)計g(r)的步數(shù)間隔。,分子模擬中徑向分布函數(shù)g(r)的算法：,第二步：統(tǒng)計g(r):,分子模擬中徑向分布函數(shù)g(r)的算法：,第三步：系綜統(tǒng)計平均計算g(r):,2. 巨正則系綜蒙特卡羅模擬方法 (Grand Canonical MC Simulation),恒定 V, T, 和 ?，體系的粒子數(shù)發(fā)生波動；可用于預測 EOS-type的性質(zhì)，但主要是用來

46、模擬吸附過程；系統(tǒng)的微觀態(tài)的分布與正則系綜類似；建立的隨機過程須增添粒子數(shù)的隨機增減；模擬是一個等化學勢面上的Markov過程。,巨正則系綜配分函數(shù),對于原子系統(tǒng)，位型（構(gòu)型）的配分函數(shù)：,其中， s 為標度坐標，r = V1/3 。概率密度為：,Metropolis GCMC algorithm,產(chǎn)生巨正則系綜的馬爾可夫鏈的過程涉及到三種不同的隨機移動： 隨機移動一個粒子，算法與正則系綜相同； 隨機添加一個粒子 隨機刪除

47、一個粒子,三種隨機移動方式的接受概率：,采用逸度（Fugacity）：,,Panagiotopoulos in 1987 introduced a clever way to simulate coexisting phases without an interface,Gibbs 系綜 MC (GEMC),,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,

48、,,,,,Two simulation volumes,,Thermodynamic contact without physical contact,MC simulation includes moves that couple the two simulation volumes,,,,,,,,,,Particle exchange equilibrates chemical potential,Volume exchange e

49、quilibrates pressure,the coupled moves enforce mass and volume balance,尤其適用于研究純流體或混合物的相平衡問題；此方法不能用于涉及到非常稠密流體的相平衡問題；此方法能同時獲得共存相的各自密度及其組成；此方法避免了共存相界面的問題。臨界點附近由于波動很大，很難模擬，所以臨界點常采用其它方式確定。,GEMC模擬方法的特點:,GEMC 的配分函數(shù),對于原子系統(tǒng)，位

50、型（構(gòu)型）的配分函數(shù):,N=N1+N2V=V1+V2 constant T,GEMC 模擬算法：,隨機選擇一個粒子進行移動 (NVT).改變每個模擬盒子的體積，但總體積保持不變 (NTP).盒子間交換粒子 (?VT)。,Phase Behavior of Alkanes,Panagiotopoulos group,Alkane Mixtures,Siepmann group:乙烷＋庚烷Ethane + n-heptane