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

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

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

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

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

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

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

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

8、邊界的表面效應(yīng)；保證中心盒子中的粒子數(shù)恒定。只需要跟蹤中心盒子中各粒子的運(yùn)動(dòng)。,周期性邊界條件的算法：,采用數(shù)學(xué)函數(shù)：,FLOOR(r/L): 返回不超過(guò)r/L的最大整數(shù),FLOOR (4.8) has the value 4.FLOOR (-5.6) has the value -6.,采用數(shù)學(xué)函數(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),定義：,中心元胞中的一個(gè)粒子只與此元胞中的其它N－1個(gè)粒子，或它們的最近鄰影像發(fā)生相互作用。,適用條件：,粒子間相互作用勢(shì)能的截?cái)嗑嚯x必須不大于模擬中心元胞長(zhǎng)度的一半。,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,此兩粒子與中心粒子的距離相等，但是：黑色球發(fā)生作用綠色球不發(fā)生作用,,此兩粒子是與中心

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

11、粒子間的相互作用主要來(lái)自于截?cái)喾秶鷥?nèi)，而范圍外的貢獻(xiàn)很小，可忽略不計(jì)。,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,截?cái)喾秶鷥?nèi)的相互作用,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,截?cái)鄤?shì)能函數(shù)的形式：,簡(jiǎn)單截?cái)鄤?shì)能函數(shù)(Truncated Potential)：,缺點(diǎn)：,rc: 截?cái)嗑嚯x或半徑,勢(shì)能在截?cái)嗵幉贿B續(xù),當(dāng)一對(duì)分子穿越邊界時(shí)，總能量不守恒。 分子間力在截?cái)嗵帪闊o(wú)窮大，MD運(yùn)

12、動(dòng)過(guò)程不穩(wěn)定。,,忽略截?cái)喟霃街獾乃凶饔?位移截?cái)鄤?shì)能函數(shù)(Shifted and Truncated Potential)：,,缺點(diǎn)：,分子間力仍然在截?cái)嗵幉贿B續(xù)。,優(yōu)點(diǎn)：,勢(shì)能在截?cái)嗵庍B續(xù)，但不影響分子間力的大小 分子間力在截?cái)嗵幉粸闊o(wú)窮大,截?cái)鄤?shì)能函數(shù)的形式：,常用于MC和MD模擬中,,,,位移－力截?cái)鄤?shì)能函數(shù)(Shifted-Force Potential)：,常用于MD模擬中,優(yōu)點(diǎn)：,勢(shì)能和分子間力均在截?cái)嗵庍B續(xù),截?cái)鄤?shì)

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

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

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

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

17、te Carlo模擬的基礎(chǔ)，服從任意分布的隨機(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ì)算機(jī)內(nèi)存采用物理方法，如利用電子線路的熱噪聲等； 昂貴而且不便重復(fù),:,?偽隨機(jī)數(shù)（Pseudo-Random Number),遞推到一定次數(shù)后，出現(xiàn)周期性的重復(fù)現(xiàn)象。,利用數(shù)學(xué)遞推

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

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

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

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

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

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

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

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

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

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

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

29、總維里項(xiàng)(Virial)；,Plrc：截?cái)喟霃酵鈱?duì)壓力的長(zhǎng)程校正,,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ùn)行作業(yè)題目READ (*,*) NStep ! 運(yùn)行步數(shù)READ (*,*) Iprint ! 打印步數(shù)READ (*,*) Isave ! 保存步數(shù)READ (*,*) Iratio ! 調(diào)整步數(shù)READ (*,‘(A)’) CNFile ! 位型文件READ (*,*) Dens

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

33、 ！判斷粒子發(fā)生重疊時(shí)的距離Rcut = Rcut * Sigma ！截?cái)喟霃紻Rmax = 0.15 * Sigma ！隨機(jī)移動(dòng)的最大距離DensLJ = DensDens = Dens / ( Sigma ** 3 )IF ( Rcut .GT. 0.5 ) Stop ' Cut-Off Too Large ',,模擬盒子的邊長(zhǎng)為1,,,

34、** Read Initial Configuration *Call ReadCN（CNFile）,正則系綜MC模擬程序F11講解（LJ, NVT）：,初始位型：,參閱程序F23,Call FCC,需要自己給定所有粒子初始位置,面心立方 (face-centered cubic, FCC)：,正則系綜MC模擬程序講解（LJ, NVT）：,長(zhǎng)程校正：,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算法,算法：勢(shì)能的計(jì)算,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 ),隨機(jī)移動(dòng),** 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,算法：步長(zhǎng)的調(diào)整,分子模擬中徑向分布函數(shù)g(r)的計(jì)算：,計(jì)算表達(dá)式：,物理意義：流體中距一個(gè)分子為 r 處出現(xiàn)另一個(gè)分子的幾率密度，它反映流體中短程有序的特點(diǎn)。,,,,,r,,r+?r,,,在模擬中通常在盒子長(zhǎng)度一半的范圍內(nèi)，考察g(r)隨距離的變化。,分子模擬中徑向分布函數(shù)g(r)的算法：,第一步：計(jì)算球殼層間的距離，并初始化一些變量：,Ngr: g(r)的統(tǒng)計(jì)次數(shù),Delg: 球殼層間的距離,nhi

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

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

47、一個(gè)粒子,三種隨機(jī)移動(dòng)方式的接受概率：,采用逸度（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,尤其適用于研究純流體或混合物的相平衡問(wèn)題；此方法不能用于涉及到非常稠密流體的相平衡問(wèn)題；此方法能同時(shí)獲得共存相的各自密度及其組成；此方法避免了共存相界面的問(wèn)題。臨界點(diǎn)附近由于波動(dòng)很大，很難模擬，所以臨界點(diǎn)常采用其它方式確定。,GEMC模擬方法的特點(diǎn):,GEMC 的配分函數(shù),對(duì)于原子系統(tǒng)，位

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