分子模擬教程.ppt

1、第三章 分子模擬方法,1. 蒙特卡羅(Monte Carlo)方法基礎(chǔ) 2. 分子動力學(xué)(Molecular Dynamics)方法基礎(chǔ) 3. 程序講解,本章具體講解內(nèi)容,掌握分子模擬方法的必備知識：,編程技能 (Fortran or C/C+) 統(tǒng)計物理學(xué)（統(tǒng)計力學(xué)）： 統(tǒng)計物理學(xué)基礎(chǔ)； 系綜原理； 非平衡統(tǒng)計力學(xué)基礎(chǔ)； 漲落理論 分子熱力學(xué) ：分子間相互作用理論； 分布函數(shù)理論 氣體分子運(yùn)動論 其它,分子模擬的目的:,為什么要進(jìn)行分子模擬？,將分子聚集體的性質(zhì)與如下方面相聯(lián)系： 分子的微觀相互作用 分子聚集體的結(jié)構(gòu) 分子的動力學(xué)過程,分子模擬對實(shí)驗(yàn)進(jìn)行補(bǔ)充，使我們能夠： 預(yù)測現(xiàn)有或新材料的

2、性質(zhì) 在分子水平研究宏觀現(xiàn)象 獲得實(shí)驗(yàn)無法或難以發(fā)現(xiàn)的東西,什么是計算機(jī)分子模擬方法？,分子模擬的定義： 統(tǒng)計力學(xué)基本原理出發(fā)，將一定數(shù)量的分子輸入計算機(jī)內(nèi)進(jìn)行分子微觀結(jié)構(gòu)的測定和宏觀性質(zhì)的計算。,按照獲得微觀態(tài)的方法不同，分子模擬分為： 蒙特卡羅方法 (Monte Carlo, MC) 分子動力學(xué)方法 (Molecular Dynamics, MD) (3) 混合方法 (hybrid method，HM),計算機(jī)分子模擬的發(fā)展歷史：,1. 蒙特卡羅方法（MC）1953 Metropolis, Ulam, Rosenbluth and Tell， Los Alamos National Lab

3、 Monte Carlo simulation of hard sphere. 2. 分子動力學(xué)方法（MD）1957Alder and Wainwrigth， Livermore Lab Molecular dynamics simulation of hard spheres.,微觀與宏觀,分子模擬在微觀尺度與實(shí)驗(yàn)室的宏觀世界之間起著橋梁的作用：,給定分子間的相互作用 “準(zhǔn)確”預(yù)測研究體系的性質(zhì),MC與MD的區(qū)別：,MC: 構(gòu)型平均，不包含動力學(xué)部分； 利用概率行走產(chǎn)生微觀態(tài)。 MD: 時間平均，產(chǎn)生動力學(xué)性質(zhì)； 利用運(yùn)動軌線隨時間的變化來產(chǎn)生一系列微觀態(tài)。,計算機(jī)分子模擬的發(fā)展歷史（續(xù)）：

4、,從上個世紀(jì)九十年代初期以來，計算機(jī)模擬技術(shù)得到了飛速發(fā)展，主要基于三個方面的發(fā)展： 分子力場的發(fā)展（基石） （Amber，OPLS、Compass） 原子間的鍵長、鍵角、分子間的內(nèi)聚能等 模擬算法（途徑） 計算機(jī)硬件（工具）,HPCx,計算機(jī)分子模擬的特點(diǎn)：,原子水平的模擬 計算機(jī)實(shí)驗(yàn) 檢驗(yàn)理論、篩選實(shí)驗(yàn) 科學(xué)研究中的第三種方法,分子模擬中涉及的幾個基本概念：,模擬計算盒子或模擬胞腔,Simulation box (cell),裝有一定數(shù)目流體分子的研究對象，它是我們要研究的宏觀體系的縮微模型。,立方形胞腔,周期性邊界條件(Periodic boundary condition, PBC),

5、在小體系中，邊界效應(yīng)總是很顯著。 在包含1000個原子的簡單立方晶體中488個原子處于邊界上。 在包含1000000個原子的簡單立方晶體中仍然有 6%的原子在邊界上。,在模擬中，考慮具有真實(shí)邊界的對象，不切合實(shí)際： 增強(qiáng)了有限尺寸效應(yīng) 人為造成的邊界會影響流體的性質(zhì),當(dāng)某個粒子運(yùn)動出模擬盒子的某一邊界時，另外一個影像粒子從另一對立邊界進(jìn)入到此盒子中。,周期性邊界條件(Periodic boundary condition, PBC),本體體系的近似：中心盒子在X，Y和Z方向無限擴(kuò)展； 消除人為形成邊界的表面效應(yīng)； 保證中心盒子中的粒子數(shù)恒定。 只需要跟蹤中心盒子中各粒子的運(yùn)動。,周期性邊界條件

6、的算法：,采用數(shù)學(xué)函數(shù)：,FLOOR(r/L): 返回不超過r/L的最大整數(shù),FLOOR (4.8) has the value 4. FLOOR (-5.6) has the value -6.,采用數(shù)學(xué)函數(shù)：,r/L0, ANINT(r/L) = AINT(r/L+0.5),r/L0, ANINT(r/L) = AINT(r/L-0.5),周期性邊界條件的算法：,y,最小影像轉(zhuǎn)化原理(Minimum image convention),定義：,中心元胞中的一個粒子只與此元胞中的其它N1個粒子，或它們的最近鄰影像發(fā)生相互作用。,適用條件：,粒子間相互作用勢能的截斷距離必須不大于模擬中心元胞長

7、度的一半。,此兩粒子與中心粒子的距離相等，但是： 黑色球發(fā)生作用綠色球不發(fā)生作用,此兩粒子是與中心原子相互作用的最近鄰影像,最小影像轉(zhuǎn)化原理的算法：,采用數(shù)學(xué)函數(shù)：,r/L0, ANINT(r/L) = AINT(r/L+0.5),r/L0, ANINT(r/L) = AINT(r/L-0.5),截斷勢能(Truncating the Potential),本體體系采用周期性邊界條件描述： 不可能將所有粒子與它們影像粒子間的相互作用全都計算。 必須在不大于中心盒子長度的一半處進(jìn)行截斷，以便與最小影像轉(zhuǎn)化原理一致。 粒子間的相互作用主要來自于截斷范圍內(nèi)，而范圍外的貢獻(xiàn)很小，可忽略不計。,截斷范圍

8、內(nèi)的相互作用,截斷勢能函數(shù)的形式：,簡單截斷勢能函數(shù)(Truncated Potential)：,缺點(diǎn)：,rc: 截斷距離或半徑,勢能在截斷處不連續(xù),當(dāng)一對分子穿越邊界時，總能量不守恒。 分子間力在截斷處為無窮大，MD運(yùn)動過程不穩(wěn)定。,忽略截斷半徑之外的所有作用,位移截斷勢能函數(shù)(Shifted and Truncated Potential)：,缺點(diǎn)：,分子間力仍然在截斷處不連續(xù)。,優(yōu)點(diǎn)：,勢能在截斷處連續(xù)，但不影響分子間力的大小 分子間力在截斷處不為無窮大,截斷勢能函數(shù)的形式：,常用于MC和MD模擬中,位移力截斷勢能函數(shù)(Shifted-Force Potential)：,常用于MD模擬中

9、,優(yōu)點(diǎn)：,勢能和分子間力均在截斷處連續(xù),截斷勢能函數(shù)的形式：,截斷勢能函數(shù)的對比：,位移力截斷勢能,簡單截斷勢能函數(shù),一、Monte Carlo模擬方法基礎(chǔ)：,亦稱統(tǒng)計模擬或隨機(jī)抽樣方法，statistical simulation method 利用隨機(jī)數(shù)進(jìn)行數(shù)值模擬的方法,Monte Carlo名字的由來：,是由Metropolis在二次世界大戰(zhàn)期間提出的：Manhattan計劃，研究與原子彈有關(guān)的中子輸運(yùn)過程；,Nicholas Metropolis (1915-1999),Monte-Carlo, Monaco,投硬幣，擲骰子,Monte Carlo方法計算Pi值,隨機(jī)數(shù)的定義和特性,

10、什么是隨機(jī)數(shù)？,單個的數(shù)字不是隨機(jī)數(shù);,是指一個數(shù)列，其中的每一個體稱為隨機(jī)數(shù)，其值與數(shù)列中的其它數(shù)無關(guān)；,在一個均勻分布的隨機(jī)數(shù)中，每一個體出現(xiàn)的概率是均等的；,例如：在0,1區(qū)間上均勻分布的隨機(jī)數(shù)序列中，0.00001與0.5出現(xiàn)的機(jī)會均等,隨機(jī)數(shù)應(yīng)具有的基本特性,隨機(jī)數(shù)序列應(yīng)是獨(dú)立的、互不相關(guān)的(uncorrelated)：,即序列中的任一子序列應(yīng)與其它的子序列無關(guān)；,長的周期(long period)：,均勻分布的隨機(jī)數(shù)應(yīng)滿足均勻性(Uniformity)：,隨機(jī)數(shù)序列應(yīng)是均勻的、無偏的，即：如果兩個子區(qū)間的“面積”相等，則落于這兩個子區(qū)間內(nèi)的隨機(jī)數(shù)的個數(shù)應(yīng)相等。,例如：對0,1)區(qū)間

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

12、，服從任意分布的隨機(jī)數(shù)序列可以用0,1區(qū)間均勻分布的隨機(jī)數(shù)序列作適當(dāng)?shù)淖儞Q或舍選后求得。,(0,1) 2 -1(-1,1),利用隨機(jī)數(shù)表，如Tippett于1972年發(fā)表的隨機(jī)數(shù)表； 占用太多的計算機(jī)內(nèi)存 采用物理方法，如利用電子線路的熱噪聲等； 昂貴而且不便重復(fù),:,偽隨機(jī)數(shù)（Pseudo-Random Number),遞推到一定次數(shù)后，出現(xiàn)周期性的重復(fù)現(xiàn)象。,利用數(shù)學(xué)遞推公式,一旦公式和初值定下來，整個隨機(jī)數(shù)序列便被確定下來，而且每一個隨機(jī)數(shù)只被它前面的那個數(shù)唯一確定，因此這類隨機(jī)數(shù)并不是真正的隨機(jī)數(shù)。,Monte Carlo方法基本思想,當(dāng)所求的問題是某種事件出現(xiàn)的概率，或是某個隨機(jī)變量

13、的期望值時，它們可以通過某種“隨機(jī)試驗(yàn)”的方法，得到這種事件出現(xiàn)的頻率和概率，或者得到這個隨機(jī)變量的統(tǒng)計平均值，并用它們作為問題的解。,Monte Carlo方法解決的問題,問題本身是確定性問題，要求我們?nèi)ふ乙粋€隨機(jī)過程，使該隨機(jī)過程的統(tǒng)計平均就是所求問題的解。,問題本身就是隨機(jī)過程，我們可以根據(jù)問題本身的實(shí)際物理過程來進(jìn)行計算機(jī)模擬和跟蹤，并采用統(tǒng)計方法求得問題的解。,Monte Carlo方法的特點(diǎn),計算的收斂性和收斂速度均與問題的維數(shù)無關(guān)，適合解決高維問題。,對問題的適應(yīng)能力強(qiáng)。,收斂速度僅為樣本數(shù)的-1/2次，因而計算耗時大。,Monte Carlo方法的應(yīng)用舉例：,計算積分：,常用

14、的積分方法求解：,將積分區(qū)域a,b均勻地劃分成N各分區(qū)間，則積分結(jié)果可近似地表示成：,x = (b-a)/N,簡單的Monte Carlo積分方法求解：,利用均勻分布的隨機(jī)數(shù)發(fā)生器，從a,b區(qū)間產(chǎn)生一系列隨機(jī)數(shù)xi，i=1, 2, ., N,其中 X為均勻分布，并且 Xa,b,近似求解Eg(X):,近似求解積分:,隨機(jī)抽樣,當(dāng)我們用簡單Monte Carlo計算積分時，若該函數(shù)為常數(shù)函數(shù)，g(x)=constant，則取樣數(shù)不管多少，準(zhǔn)確度為100。 如果在積分區(qū)間內(nèi)，g(x)為一平滑函數(shù)，則簡單Monte Carlo方法較為準(zhǔn)確，反之，如果g(x)的變動很劇烈，則簡單Monte Carlo方

15、法的誤差會變大。,說明：,重要性Monte Carlo抽樣方法,在 g(x) 變化劇烈時，如果以Monte Carlo方法取樣，最好依據(jù)g(x)的大小來決定取樣率。 當(dāng)|g(x)|的值較大時，對g(x)dx的貢獻(xiàn)也較大，如果沒被選中，則結(jié)果的誤差極大。 解決方式：改變 x被選中的機(jī)率，讓|g(x)| 值較大的點(diǎn)被選中的機(jī)率增加。 采用權(quán)重分布函數(shù)(Weight distribution function) w(x) :決定每個x被選中的機(jī)率。,重要性抽樣的定義：根據(jù)一定的分布形式進(jìn)行的隨機(jī)抽樣。,w(x)必須歸一化，即在積分區(qū)間內(nèi)w(x)dx=1。 由于 x 的選取已被 w(x) 扭曲，所以計

16、算積分時要把這部分還回去：若一共取樣了N個x，則積分值為：,重要性Monte Carlo抽樣方法,Metropolis Monte Carlo方法,我們所模擬的系統(tǒng)最終要達(dá)到的平衡分布是Boltzman分布：,Boltzmann 概率分布函數(shù)：,我們?nèi)绻軌虍a(chǎn)生這種分布，我們就能夠計算系統(tǒng)的大多數(shù)性質(zhì)，但這是不可能的，因?yàn)槲覀儾恢繸的值，但是對于任意兩個狀態(tài)，我們有：,可以在相空間中構(gòu)造一個馬爾科夫鏈，使相空間中的樣本點(diǎn)隨著鏈的增長逐步趨近于Boltzman分布。,一個序列x0, x1, x2, ,xn,如果對任何n都有：,則此序列是一個Markov鏈。,要求：,任何一次實(shí)驗(yàn)的結(jié)果依賴于前一

17、次的試驗(yàn)，并且近依賴于前一次的試驗(yàn)。,馬爾科夫(Markov)鏈,可以證明：通過構(gòu)造Markov鏈，體系中最終的平衡分布就是Boltzman分布,Metropolis平衡條件 (Detailed balance condition)：,平衡條件：,系統(tǒng)處于狀態(tài)X的概率正比于其Boltzman因子：,如果是對稱的:,Metropolis Monte Carlo方法的算法：,給出一個初始狀態(tài)，并計算系統(tǒng)的能量Eold,隨機(jī)產(chǎn)生一個新狀態(tài)，并計算新系統(tǒng)的能量Enew,如果E(EnewEold)0, 則接受新狀態(tài)并回到 b),如果E(EnewEold)0, 則計算Boltzman因子：,在(0, 1)

18、區(qū)間上產(chǎn)生一個均勻分布的隨機(jī)數(shù) ；,如果,則接受新狀態(tài)并回到b),否則保留原值并回到b),1. 正則系綜蒙特卡羅模擬方法(Canonical MC Simulation),具有確定的粒子數(shù)N、溫度T和體積V,對于含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,隨機(jī)產(chǎn)生一個新狀態(tài)，并計算新系統(tǒng)的能量Unew,如果U(UnewUold)0,

19、 則接受新狀態(tài)并回到 b),如果U(UnewUold)0, 則計算Boltzman因子：,在(0, 1)區(qū)間上產(chǎn)生一個均勻分布的隨機(jī)數(shù) ；,如果,則接受新狀態(tài)并回到b),否則保留原值并回到b),正則系綜MC模擬算法的組織：,正則系綜MC模擬算法的流程圖：,給定每個分子的初始位置，ri(0),隨機(jī)選取一個分子，并隨機(jī)移動到新的位置,計算移動前后的系統(tǒng)能量變化U,拒絕移動,U0 ?,Exp(U) (0,1)?,統(tǒng)計系統(tǒng)的熱力學(xué)性質(zhì)及其它物理量,統(tǒng)計性質(zhì)不變？,打印結(jié)果，結(jié)束,Yes,No,No,No,接受移動,Yes,Yes,大約循環(huán)107到108次,Monte Carlo模擬中幾個熱力學(xué)量的計算

20、：,N個粒子系統(tǒng)中的總勢能：,假設(shè)采用截斷勢能函數(shù)：,Uc：截斷范圍內(nèi)的總勢能； Ulrc：截斷半徑外對勢能的長程校正(Long-range correction),對于LJ流體：,含N個粒子系統(tǒng)中的壓力：,Wc：截斷范圍內(nèi)的總維里項(xiàng)(Virial)；,Plrc：截斷半徑外對壓力的長程校正,Read simulation parameters,Start,Initialize positions of all particles,New simulation?,Read old configuration,Monte Carlo loop,Stop,yes,no,Monte Carlo loo

21、p Subroutine,Start,Stop,Trial move,Satisfy Metropolis rule?,Accept the trial move,Update energy and virial,Sample the pressure,End of simulation?,yes,no,yes,no,Main program,正則系綜MC模擬程序基本結(jié)構(gòu)：,正則系綜MC模擬程序F11講解（LJ, NVT）：,* READ INPUT DATA *,初始狀態(tài)：,READ (*,(A) Title ! 運(yùn)行作業(yè)題目 READ (*,*) NStep ! 運(yùn)行步數(shù) READ (*,

22、*) Iprint ! 打印步數(shù) READ (*,*) Isave ! 保存步數(shù) READ (*,*) Iratio ! 調(diào)整步數(shù) READ (*,(A) CNFile ! 位型文件 READ (*,*) Dens ! 對比密度 READ (*,*) Temp ! 對比溫度 READ (*,*) Rcut ! 對比截斷半徑,無因次量：,正則系綜MC模擬程序F11講解（LJ, NVT）：,量綱變換：,Beta = 1.0 / Temp Sigma = ( Dens / Real ( N ) ) * ( 1.0 / 3.0 ) Rmin = 0.70 * Sigma ！判斷粒子發(fā)生重疊時的距離 R

23、cut = Rcut * Sigma ！截斷半徑 DRmax = 0.15 * Sigma ！隨機(jī)移動的最大距離 DensLJ = Dens Dens = Dens / ( Sigma * 3 ) IF ( Rcut .GT. 0.5 ) Stop Cut-Off Too Large ,模擬盒子的邊長為1,* Read Initial Configuration * Call ReadCN（CNFile）,正則系綜MC模擬程序F11講解（LJ, NVT）：,初始位型：,參閱程序F23,Call FCC,需要自己給定所有粒子初始位置,面心立方 (face-centered cubic, FCC)

24、：,正則系綜MC模擬程序講解（LJ, NVT）：,長程校正：,Sr3 = ( Sigma / Rcut ) * 3 Sr9 = Sr3 * 3 Vlrc12 = 8.0 * Pi * DensLJ * Real ( N ) * Sr9 / 9.0 Vlrc6 = - 8.0 * Pi * DensLJ * Real ( N ) * Sr3 / 3.0 Vlrc = Vlrc12 + Vlrc6 Wlrc12 = 4.0 * Vlrc12 Wlrc6 = 2.0 * Vlrc6 Wlrc = Wlrc12 + Wlrc6,算法：能量求和：,Call Sumup ( Rcut, Rmin, Sig

25、ma, Ovrlap, V, W ) If ( Ovrlap ) Stop Overlap In Initial Configuration Vs = ( V + Vlrc ) / Real ( N ) Ws = ( W + Wlrc ) / Real ( N ) Ps = Dens * Temp + W + Wlrc Ps = Ps * Sigma * 3,* Check For Acceptance * DeltV = Vnew - Vold DeltW = Wnew - Wold DeltVb = Beta * DeltV If ( DeltVb .Lt. 75.0 ) Then If

26、( DeltV .Le. 0.0 ) Then V = V + DeltV W = W + DeltW RX(i) = RXInew . Acatma = Acatma + 1.0 Elseif ( Exp ( - DeltVb ) .Gt. Ranf ( Dummy ) ) Then V = V + DeltV W = W + DeltW RX(i) = RXInew Acatma = Acatma + 1.0 Endif Endif Acm = Acm + 1.0,算法：Metropolis算法,算法：勢能的計算,DO 100 J = 1, N IF ( I .NE. J ) Then R

27、XIJ = RXI - RX(J) RXIJ = RXIJ - Anint ( RXIJ ) RIJSQ = RXIJ * RXIJ + RYIJ * RYIJ + RZIJ * RZIJ IF ( RIJSQ .LT. RcutSQ ) THEN SR2 = SIGSQ / RIJSQ SR6 = SR2 * SR2 * SR2 VIJ = SR6 * ( SR6 - 1.0 ) WIJ = SR6 * ( SR6 - 0.5 ) V = V + VIJ W = W + WIJ ENDIF ENDIF 100 Continue V = 4.0 * V W = 48.0 * W / 3.0,最

28、小影像原理,算法：Metropolis算法,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.0 * Ranf ( DUMMY ) - 1.0 ) * DRmax . RXInew = RXInew - Anint ( RXInew ) * C

29、alculate The Energy Of I In The New Configuration * Call Energy ( RXInew, RYInew, RZInew, I, Rcut, Sigma,Vnew, Wnew ),隨機(jī)移動,* Adjust Maximum Displacement * Ratio = Acatma / Real ( N * Iratio ) If ( Ratio .Gt. 0.5 ) Then DRmax = DRmax * 1.05 Else DRmax = DRmax * 0.95 Endif Acatma = 0.0 Endif If ( Mod

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

31、):,分子模擬中徑向分布函數(shù)g(r)的算法：,第三步：系綜統(tǒng)計平均計算g(r):,2. 巨正則系綜蒙特卡羅模擬方法 (Grand Canonical MC Simulation),恒定 V, T, 和 ，體系的粒子數(shù)發(fā)生波動； 可用于預(yù)測 EOS-type的性質(zhì)，但主要是用來模擬吸附過程； 系統(tǒng)的微觀態(tài)的分布與正則系綜類似； 建立的隨機(jī)過程須增添粒子數(shù)的隨機(jī)增減； 模擬是一個等化學(xué)勢面上的Markov過程。,巨正則系綜配分函數(shù),對于原子系統(tǒng)，位型（構(gòu)型）的配分函數(shù)：,其中， s 為標(biāo)度坐標(biāo)，r = V1/3 。概率密度為：,Metropolis GCMC algorithm,產(chǎn)生巨正則系綜的馬爾可夫鏈的過程涉及到三種不同的隨機(jī)移動： 隨機(jī)移動一個粒子，算法與正則系綜相同； 隨機(jī)添加一個粒子 隨機(jī)刪除一個粒子,三種隨機(jī)移動方式的接受概率：,采用逸度（Fugacity）：,Panagiotopoulos in 1987 introduced a clever way to simulate coexisting phases without an interface,Gibbs 系綜 MC (GEMC),Two simulation volumes,Thermodynamic contact without physical contact,MC sim