物理問題的計(jì)算機(jī)模擬方法

1、碩士研究生課程物理問題的計(jì)算機(jī)模擬方法講義適用專業(yè)：凝聚態(tài)物理、材料物理與化學(xué)、理論物理、光學(xué)工程學(xué) 時(shí)：3040學(xué)時(shí)參考教材：德著，秦克誠(chéng)譯，理論物理中的計(jì)算機(jī)模擬方法，北京大學(xué)出 版社，1996。荷Frenkel & Smit著，汪文川 等譯，分子模擬一從算法到應(yīng)用，化學(xué)工業(yè)出版社，2002。and Computer Simulation of Liquids, Clarendon Press, Oxford, 1989.Molecular Modelling： Principles and Applications, Addison Wesley Longman, England, 19

2、96.德D.羅伯 著，計(jì)算材料學(xué)，化學(xué)工業(yè)出版社，2002。6.英B. Chopard & Michel Droz著，物理系統(tǒng)的元胞自動(dòng)機(jī)模擬，祝玉學(xué)，趙學(xué)龍 譯，清華大學(xué)出版社，2003。第一章計(jì)算機(jī)模擬方法概論序言熱力學(xué)系統(tǒng)物理量的統(tǒng)計(jì)平均分子動(dòng)力學(xué)方法模擬的基本思想蒙特卡羅方法模擬的基本思想元胞自動(dòng)機(jī)模擬的基本思想簡(jiǎn)要的發(fā)展歷程簡(jiǎn)單元胞自動(dòng)機(jī)：奇偶規(guī)則元胞自動(dòng)機(jī)的一般定義第二章確定性模擬方法一分子動(dòng)力學(xué)方法（MD）分子動(dòng)力學(xué)方法微正則系綜分子動(dòng)力學(xué)方法正則系綜分子動(dòng)力學(xué)方法 等溫等壓系綜分子動(dòng)力學(xué)方法第三章 隨機(jī)性模擬方法一蒙特卡羅方法（MC）預(yù)備知識(shí)布朗動(dòng)力學(xué)（BD）蒙特卡羅方法微正則

3、系綜蒙特卡羅方法正則系綜蒙特卡羅方法等溫等壓系綜蒙特卡羅方法巨正則系綜蒙特卡羅方法第四章離散性模擬方法一原胞自動(dòng)機(jī)（CA）引言元胞自動(dòng)機(jī)模擬 *元胞自動(dòng)機(jī)模擬的應(yīng)用第一章計(jì)算機(jī)模擬方法概論序言1.什么是計(jì)算機(jī)模擬Simulation Modelling2.為什么要進(jìn)行計(jì)算機(jī)模擬3.常用的計(jì)算機(jī)模擬方法確定性模擬方法：MD模擬Molecular Dynamics隨機(jī)性模擬方法：MC模擬Monte Carlo離散性模擬方法：CA模擬Cellular Automata熱力學(xué)系統(tǒng)物理量的統(tǒng)計(jì)平均描述系統(tǒng)的坐標(biāo)(自由度)：x(t) = xi (t),x2 (t),XN (t)系統(tǒng)的物理量：A(X(t)1

4、 ,時(shí)間平均(1-1)A =L f tA(x(t)dt分子動(dòng)力學(xué)(MD)模擬I。t0系綜平均= I!期)亍(H(X)dX蒙特卡羅(MC)模擬(1-2)=J A(x)p (x)dxQp (x) = 1 f (H (x)一 分布函數(shù)(幾率密度函數(shù)) Z(1-3)Z = j f (H(x)dx配 分 函 數(shù)Q(1-4)Q一相空間H (x)系統(tǒng)的哈密頓函數(shù)對(duì)于處于平衡態(tài)的系統(tǒng)，可以證明：A = A8對(duì)于實(shí)際的有限時(shí)間內(nèi)的平均，則有A x A實(shí)際模擬的系統(tǒng)大小也是有限的：有限的粒子數(shù)N或有限的系統(tǒng)限度L對(duì)統(tǒng)計(jì)平均結(jié)果有影響。分子動(dòng)力學(xué)(MD)方法模擬的基本思想1.基本原理系統(tǒng)：N個(gè)粒子，體積V,粒子質(zhì)量

5、為m描述一個(gè)粒子運(yùn)動(dòng)狀態(tài)的自由度：(r., p.)(p_=mv_)相空間：6N維，相空間中的一點(diǎn)的坐標(biāo)Xn=rN, (mvN)rN=(r, r ,，r ),vn=(v, v,，v)12N12N粒子間的相互作用勢(shì)：U(rN)=U(r1, r2,，匚)二寸u(r )y決定系統(tǒng)相軌跡Xn (t)的運(yùn)動(dòng)方程：drdv/ =mdf f (r J(i = 1,2,.,N)Xn(0) = Xn (初始條件) 0、加上邊界條件(周期性)(1-5)物理量A的宏觀值，由A(Xn)的時(shí)間平均獲得，即如)=1 / AX(tf)dtf(離散情況：A(t) = - A)oi=1對(duì)于平衡態(tài)：A = lim A (t) t

6、s實(shí)際模擬時(shí)間總是有限的，模擬時(shí)間的長(zhǎng)短可通過判斷時(shí)間的增加對(duì)平均值的影響來確定，當(dāng)繼續(xù)增加時(shí)間帶來的平均值得變化在允許的誤差范圍之內(nèi)時(shí)，即可認(rèn)為模擬足夠長(zhǎng)了。2.計(jì)算步驟運(yùn)動(dòng)方程：空i_ = v , m空 = -V U(rn ) dtdt即md = -V U(rn) = F(1-6)或d = F/ m(1-7)數(shù)值求解：用差分近似表示微分(采用不同的差分格式，可得到不同的算法)。用顯示中心差分格式，將(7)式寫為 TOC o 1-5 h z 竺 r r (t + At) - 2r (t) + r (t-At )/(At) 2(1-8)由(7)和(8)式可得：r (t + At) = 2r (

7、t) - r (t - At) + (At)2F /m(1-9)iiii第一步：由(9)式計(jì)算第i個(gè)粒子在t+At時(shí)刻的位置坐標(biāo)。要啟動(dòng)計(jì)算，我們必須要知道最初兩點(diǎn)ri(0)和匚(七)第二步：對(duì)不同時(shí)刻t = At , 2At, 3At,，LAt(t0 =0)計(jì)算物理量A(r1(lAt), r2(lAt),，rN(lAt)(l = 1,，N)第三步：計(jì)算物理量A的平均值1= lim-1 A(r (W), r (lAt), , r (lAt)lM -12Nl =1L的大小由繼續(xù)增大L而不變(或變化在誤差范圍內(nèi))來 確定。蒙特卡羅(MC)方法模擬的基本思想1.基本原理以正則系綜(T, V, N)為

8、例正則分布：p = Z e-區(qū)正則配分函數(shù)：z = 1一I e推dQ(dQ = drNdpn)N !h 3 n系統(tǒng)能量：E = E + U(rN) = Z 虹 + U(rn) s p2mi i物理量：A(rN)= A(r1, r1,， ) 系綜平均:J A(r n)p1de=N!h3ns1J A(r n )e - es / kT dr Ndp n N!h3NZ1-S 略 / kT ,J A(rn )e-u(rn)/kTdrn J e . 2m,dpn N!h3NZ!J A(rn )e-u(rn)/kTdrnQNQ = J e -u (rN) / kT dr n（位形積分）(1-11)用MC方法

9、計(jì)算上述多維積分。2.計(jì)算步驟（1）劃分原胞N個(gè)粒子一3N個(gè)空間自由度，3N維空間劃分成s個(gè)相等的原胞（s1）注意：由于積分中不含動(dòng)量，所以我們只需要在位置空間積分，而不 需要在相空間中積分。當(dāng)系統(tǒng)的代表點(diǎn)落入第i個(gè)原胞時(shí)，則認(rèn)為系統(tǒng)處在狀態(tài)i,因此，s 為系統(tǒng)可能的微觀狀態(tài)數(shù)目。于是，積分（10）和（11）可近似表述為41 y .= A e -ui/ kTQ N ie-Ui / kT(1-13)(2)建立馬爾可夫（MapkoB）過程（鏈）將s個(gè)狀態(tài)看作一組隨機(jī)事件馬爾可夫鏈：從狀態(tài)iTj狀態(tài)j（ I）的概率pij ,只與 Zp =1, i=1, 2,，sj若j經(jīng)歷n步到達(dá).其概率表示為p（n

10、），存在極限概率lim p( n)= u ns j(j=1,2,，s)u 0, Zu. = 1ju.為系統(tǒng)處在狀態(tài)j的概率。于是，沿?zé)o限長(zhǎng)的馬爾可夫鏈，物理量A的平均值可寫為=Z Au (=ZZ Ap )i ii ij(1-14)選取-1 e -u,/kT ,N則（14）式為A的正則系綜平均值。（3）抽樣方法采用怎樣的抽樣方法所構(gòu)成的馬爾可夫鏈能得到上述平均值粒子位置坐標(biāo)：x a （i）1粒子編號(hào)： =1,2,N坐標(biāo)的三個(gè)方向:=1, 2, 3系統(tǒng)狀態(tài)：i=1,2,，s給定粒子位置坐標(biāo)的變化量（小于系統(tǒng)體積的限度）給定系統(tǒng)的初態(tài)i,隨機(jī)選定4個(gè)隨機(jī)數(shù)，其中三個(gè)（=1, 2,3）,且-11, 一

11、個(gè)表示粒子編號(hào)二1, 2,N,由此隨機(jī)確定粒子位置的變化：Xa(i) T Xa(j) = Xa(i) +如8111(確保xa(i) - Xa(j) 5 )11若U U，則再選一個(gè)隨機(jī)數(shù)4 （04 1 ）,若&4 e-（U廠U）/kT，則粒子保留在原位置，不發(fā)生i j的躍遷;若 g 4 eTUj-Ui)/kT，則發(fā)生i j的躍遷。由此進(jìn)行下去，則形成一個(gè)馬爾可夫鏈(或過程)，此鏈的長(zhǎng)度L(即粒 子行走的步數(shù)，遠(yuǎn)大于s)，由所計(jì)算的物理量的平均值= lim L i AlfL i=1 (1-15)不再隨鏈的加長(zhǎng)而改變來確定。由此得到的平均值即正則平均值。一般來 說，L 與 N, V, T 有關(guān)，比如

12、，N=32108, L=30005000。歸納起來，計(jì)算系統(tǒng)物理量的正則系綜平均值的具體步驟如下：第一步：給定系統(tǒng)的初始狀態(tài)(粒子的初始位置)匚和每一步的改變量 ,第二步：選擇四個(gè)隨機(jī)數(shù)，其中一個(gè)代表粒子的編號(hào)i ( 1 i N )； 另外三個(gè)表示粒子空間坐標(biāo)的改變x, y, z (-, =1,2, 3);第三步：計(jì)算粒子i的新位置rr= 3,礦,Z) = (x +5 ,y +8 ,z +8 )i i i ii x i y i z第四步：計(jì)算粒子在新舊兩個(gè)位置系統(tǒng)的能量之差A(yù)U = U(r ,r ,., r,., r ) - U(r ,r ,., r ,., r ) 12 i N12 i N第

13、五步：由AU的大小判斷粒子i是否從ri運(yùn)動(dòng)到r:若AU 0，則再選一個(gè)隨機(jī)數(shù)R（0 R 1）,如果R e _ au/燈，則J不變，返回第二步。第六步：計(jì)算A（r,r，,/,.,r ） 12 i N第七步：重復(fù)上述各個(gè)步驟，直到完成L步為止，最后利用公式（15）計(jì)算A的平均值。粒子間相互作用勢(shì)模型的選取最簡(jiǎn)單的兩種模型:（1）硬球模型I 0u (r )= 3(r a)(r b)（ 為硬球的直徑）LJ 勢(shì)u (r)=-r 6Br 12U(r ,r，,r ) = u(r )12 Nij 元胞自動(dòng)機(jī)(CA)模擬的基本思想元胞自動(dòng)機(jī)：時(shí)間和空間都離散、物理參量只取有限數(shù)值集的物理系統(tǒng)的 理想化模型cel

14、lular automata 或 cellular automaton CA簡(jiǎn)要的發(fā)展歷程自繁殖系統(tǒng)20世紀(jì)40年代，Von Neumann,構(gòu)造能解決非常復(fù)雜問題的計(jì)算機(jī)， 設(shè)想模仿人腦的行為尋求與生物過程無(wú)關(guān)的情況下自繁殖機(jī)理的邏輯 抽象。根據(jù)S. Ulam的建議，Von Neumann在由元胞構(gòu)成的完全離域的框架 下處理這個(gè)問題，構(gòu)造了一個(gè)完全離散的動(dòng)力學(xué)系統(tǒng)元胞自動(dòng)機(jī)。第一個(gè)自復(fù)制元胞自動(dòng)機(jī)由二維方形網(wǎng)絡(luò)組成，有數(shù)千個(gè)基本元胞構(gòu)成的自繁殖結(jié)構(gòu)。(1) 一般機(jī)器只能構(gòu)造比自己簡(jiǎn)單的客體，而采用自復(fù)制元胞機(jī)，可 獲得一種能產(chǎn)生新的、具有同樣復(fù)雜性和功能的“機(jī)器”；(2) Von Neuma

15、nn的元胞自動(dòng)機(jī)規(guī)則具有所謂通用計(jì)算的性質(zhì)，這意 味著，存在一種元胞自動(dòng)機(jī)的初始構(gòu)型，該元胞自動(dòng)機(jī)能產(chǎn)生任何計(jì)算機(jī) 算法的解。通用計(jì)算的性質(zhì)指：用元胞自動(dòng)機(jī)演化規(guī)則能夠模擬任何計(jì)算 機(jī)流程（邏輯選擇器開關(guān)）。生命游戲機(jī)1970年，數(shù)學(xué)家John Conway生命游戲機(jī)的概念，尋找能導(dǎo)致復(fù)雜行 為的簡(jiǎn)單規(guī)則。設(shè)想一個(gè)類似于棋盤的二維方形網(wǎng)格，每個(gè)元胞可能的狀態(tài)是活（狀態(tài) 1）或死（狀態(tài)0），其更新規(guī)則是：有三個(gè)活元胞包圍的一個(gè)死元胞恢復(fù)為 活元胞；由兩個(gè)以下或三個(gè)以上活元胞包圍的活元胞因孤立或擁擠而死亡。結(jié)果表明，生命游戲機(jī)有出乎意料的豐富行為，從原“湯”中顯示出來 的復(fù)雜結(jié)構(gòu)，演變發(fā)展成為某些

16、特殊的技藝，例如，可能形成所謂的滑翔 機(jī)緊鄰元胞的特殊排列，這些元胞具有沿直線彈道穿越空間運(yùn)動(dòng)的特 性。生命游戲機(jī)也是具有計(jì)算通用性的元胞自動(dòng)機(jī)。模擬物理系統(tǒng)（1）20世紀(jì)70年代，Hardy, Pomeau和Pazzis建立了所謂的HPP 格子氣體模型，用以在質(zhì)量和動(dòng)量守恒的情況下在方形網(wǎng)格上模擬粒子的 碰撞行為。(2) 1986 年，F(xiàn)risch, Hasslacher 和 Pomeau 提出了著名的 FHP 模 型，這是在六邊形網(wǎng)格上模擬二維流體動(dòng)力學(xué)的第一個(gè)嚴(yán)格模型全離 散計(jì)算機(jī)模型替代風(fēng)洞試驗(yàn)。HPP和FHP通常稱之為格子氣自動(dòng)機(jī)（LGALattice Gas Automata）（3

17、）Ising自旋動(dòng)力學(xué)模型，20世紀(jì)80年代末，Vichniac提出Q2R 規(guī)則。（4）格子Boltzmann方法與多粒子模型格子Boltzmann方法或模型（LBM）：網(wǎng)格上定義的物理模型，在這個(gè) 網(wǎng)格上與每個(gè)格位相關(guān)聯(lián)的變量是平均粒子數(shù)，或具有一定速度粒子出現(xiàn) 的概率。該模型可以用均化或因子分解方法由元胞自動(dòng)機(jī)動(dòng)力學(xué)推導(dǎo)出來， 或者自定義，而與特定的實(shí)現(xiàn)無(wú)關(guān)。格子Boltzmann模型保持了元胞自動(dòng)機(jī)方法的微觀水平解釋，但忽略 了多體的相關(guān)函數(shù)，但這種方法已經(jīng)成為目前模擬物理系統(tǒng)中最有前途的 方法之一。在嚴(yán)格的元胞自動(dòng)機(jī)方法與較靈活的格子Boltzmann模型之間，有一 種目前處于發(fā)展中的

18、模型多粒子模型。這種模型保留量化狀態(tài)的概念， 但接受無(wú)限數(shù)值集，因此既保證了數(shù)值穩(wěn)定性（與LBM相反），又考慮了多 體相關(guān)性。在模擬物理系統(tǒng)時(shí)，大量的可能狀態(tài)提供更多的靈活性，并產(chǎn) 生小的統(tǒng)計(jì)噪聲。但多粒子動(dòng)力學(xué)更難設(shè)計(jì)，且在數(shù)值計(jì)算上比格子 Boltzmann方法更慢。簡(jiǎn)單元胞自動(dòng)機(jī)：奇偶規(guī)則簡(jiǎn)單元胞自動(dòng)機(jī)的演化規(guī)則（20世紀(jì)70年代，Edward Fredkin提出， 定義在二維方形網(wǎng)上）：網(wǎng)格的每一個(gè)格位是一個(gè)元胞，以其位置r =（i, j）來標(biāo)記，其中i 和為行和列的標(biāo)號(hào)。函數(shù)寸（r）描述每個(gè)元胞在時(shí)間t的狀態(tài)，其值為0 或1。從時(shí)間t = 0的初始條件及網(wǎng)格上給定的構(gòu)形值寸o（r）開

19、始，t = 1時(shí) 狀態(tài)按下列步驟求得：（1）對(duì)于每個(gè)格位r，都計(jì)算出其位于東、南、西、北4個(gè)最近鄰格 位r，的wo（r）值之和。應(yīng)使系統(tǒng)在i和j兩個(gè)方向循環(huán)（如同在環(huán)面上）， 從而確定出所有格位的W （r，）計(jì)算值。（2）如果這個(gè)和值為偶數(shù)，則新狀態(tài)寸1（r）為0（白色），否則為1 （黑 色）。重復(fù)上述步驟，得到t = 2, 3, 4,的狀態(tài)。這個(gè)元胞自動(dòng)機(jī)的奇偶規(guī)則可表示為：W t 1 （i, j） = W t （i +1, j） W t （i -1, j） W t （i,j +1） W t （i, j -1）式中，符號(hào) 代表“異或”邏輯運(yùn)算，也即模2和：1 1=0 0=0，1 0=0 1=

20、1。反復(fù)迭代這個(gè)規(guī)則時(shí)，可得到非常精美的幾何圖形。元胞自動(dòng)機(jī)的一般定義定義:規(guī)整的元胞網(wǎng)格覆蓋d維空間的一部分；歸屬于網(wǎng)格的每個(gè)格位r的一組布爾變量(r, t) = 1(r,t),2(r, t),，m(r, t)給出每個(gè)元胞在時(shí)間t = 0,1, 2,的局部狀態(tài)；演化規(guī)則R= R1, R2,，RJ按下列方式指定狀態(tài)(r, t) 的時(shí)間演化過程：j(r, t+1) = R (r, t),(r + 1, t),(r + 2,t),，(r +, t)q式中r +指定從屬于元胞r的給定鄰居。q鄰居：二維元胞自動(dòng)機(jī)的兩種鄰居：(1)Von Neumann鄰居，有一個(gè)中心 元胞(要演化的元胞)和四個(gè)位于其

21、近鄰東西南北方位的元胞組成；(2) Moore鄰居，除了前面涉及的最近鄰元胞外，還包括次近鄰的4個(gè)元胞， 共九個(gè)元胞。還有一種有用的鄰居稱為Margolus鄰居，將空間劃分成2 2元胞 的鄰接單元塊，這個(gè)規(guī)則對(duì)位于單元塊內(nèi)的位置即左上、右上和左下、 右下很敏感。三種鄰居如下圖所示。邊界條件：周期性邊界條件；固定邊界條件；絕熱邊界條件；映射邊界條件。備注：確定型元胞自動(dòng)機(jī)：演化規(guī)則確定，給定的初始狀態(tài)將始終演化出同樣的式樣。概率型自動(dòng)元胞機(jī)：演化規(guī)則包含一定的隨機(jī)性，給定的初始狀態(tài)可能演化出不同的式樣。第二章 確定性模擬方法一分子動(dòng)力學(xué)方法（MD） 分子動(dòng)力學(xué)方法用分子動(dòng)力學(xué)方法模擬氣體、液體或

22、固體系統(tǒng)分子動(dòng)力學(xué)的描述形式：哈密頓描述拉格朗日描述牛頓運(yùn)動(dòng)方程描述劃分MD元胞選取適當(dāng)大小的V=L3,太大耗費(fèi)計(jì)算時(shí)間，太小不能準(zhǔn)確反映 系統(tǒng)的性質(zhì)。目的：維持系統(tǒng)恒定的密度。對(duì)平衡態(tài)系統(tǒng)，液體和氣體原胞的形狀無(wú)關(guān)緊要，但對(duì)晶體則 不然。周期性邊界條件A(r) = A(r + nL)目的：減少表面效應(yīng)相互作用勢(shì)、最小影像約定和切斷距離U(r ,r ,.,r ) = Zu(r ) + ZZ u(t r + nLl)n = (n , n , n )12 Nij i j 1123i r的相互作用忽略 不計(jì)，代價(jià)是忽略了背景。積分格式F = U(rn) = ZV u(r ) = ZF(r )iii

23、ijij(i ) j(i F(r )(i) j(1)遞推公式在第一章中求解運(yùn)動(dòng)方程時(shí)，我們是直接求解的關(guān)于粒子位置坐標(biāo)的 二階微分方程，得到的遞推公式需要知道最初的兩點(diǎn)位置才能啟動(dòng)計(jì)算， 但在實(shí)際計(jì)算中，我們常常是給出最初的位置和速度，于是，我們可通過 選取一定的差分格式，有運(yùn)動(dòng)方程(16)得到關(guān)于位置和速度的遞推公式。|rn)-rn+i)或；r(t)tr(t+h)已知(0)=roI V(n) T V(n+1)v(t) t v(t + h)v(0) = vk ii0遞推公式的具體形式取決于差分各式的選取。時(shí)間步長(zhǎng)選擇時(shí)間步長(zhǎng)的原則：在保證計(jì)算精度的前提下，盡量節(jié)省計(jì)算時(shí)間。實(shí)例：Ar原子系統(tǒng)，

24、LJ勢(shì)時(shí)間步長(zhǎng)取為h 10-2(約化時(shí)間步長(zhǎng))10- 14S (=100ps=10fs)約束條件保證能量、動(dòng)量或角動(dòng)量守恒。減小計(jì)算誤差的技巧數(shù)值計(jì)算不可避免有誤差，與誤差有關(guān)的因素主要有:差分格式 時(shí)間步長(zhǎng)切斷半徑最小影像約定等等計(jì)算熱力學(xué)量A = lim(t t )j A(rn (t), v(t),V(t)dt t 8t0若系統(tǒng)的NVE恒定，則A = (微正則系綜平均一 microcannonical ensemble average)若系統(tǒng)的NVT恒定，則A =(正則 系綜平均一cannonical ensemble average)在實(shí)際計(jì)算中，往往還需要要求系統(tǒng)的總動(dòng)量P=0或恒定，

25、因整 個(gè)系統(tǒng)不受外界的作用。溫度與能量均分定理-1_、E = kT(3N N )k 2cN一總粒子數(shù)，N一約束條件個(gè)數(shù)，一般情況下，N NP=0Nc = 3位勢(shì)切斷誤差與修正g(r) 一 對(duì)關(guān)聯(lián)函數(shù)(pair correlation function)p =N 一粒子數(shù)密度Vpg(r)dr一當(dāng)原點(diǎn)位置上一個(gè)粒子時(shí)，在r附近dr內(nèi)發(fā)現(xiàn)有一個(gè)粒子的幾率。對(duì)氣體或液體，p (r) = p (r)(各向同性)平均每個(gè)粒子的能量 =4 J u (r) pg (r )dr =上 J u (r) pg (r )r 2 drdQN 22=2 p兀 J u (r) g (r )r 2 dr = 2 p兀 J 5

26、 (r) g (r )r 2 dr + 2 p兀u (r) g (r )r 2 dr00rc機(jī) 2p兀 Jrc u(r)g(r)r2drN 0切斷誤差修正為:U = 2p兀 Mu(r)g(r)r2drrc6.計(jì)算機(jī)模擬的組織初始化給定粒子的初始位置、初始速度等;趨衡從初始轉(zhuǎn)臺(tái)開始，按運(yùn)動(dòng)方程要求的規(guī)律，從非平衡態(tài)趨向平衡態(tài)；投產(chǎn)計(jì)算物理量的統(tǒng)計(jì)平均值。 微正則系綜分子動(dòng)力學(xué)方法微正則系綜：NVE恒定，且P = 0 (分子動(dòng)力學(xué)方法的自然選擇)1. Verlet 算法運(yùn)動(dòng)方程的顯示中心差分格式r (t + kt) = 2r (t) - r (t -At) + (At )2 F / m iiiiv

27、. (t) = (r (t + At) - r (t - At) / 2At由可得到如下遞推公式(2-2)rn+i = 2rn rn-i + Fnh2 /m vn = (rn+1 rni) / 2h此即Verlet算法(A2),其特點(diǎn)是:要啟動(dòng)此算法，必須已知兩點(diǎn)r。, n，所以又稱為二步法;rn和vn不能同步算出，第n+1步算出的第n步的速度。2. Verlet算法的速度形式r (t + At)機(jī) r (t) + 空也 At +1 竺也(At)211 dt 2 dt 21 F A=r (t) + v. (t)At + i-(At)2F(t)v (t + At) = v (t) + 里 At_

28、m.、, A、, 、 F(t + At) Av (t + At) = v(t) +At、m(向前差分)n v(t + At) = v(t) +(向后差分)2m由此可得如下遞推公式:r n+1 = r n + hvn + h 2 F n /2mVn+1 = Vn + h( F n+1 + Fn )/2mi iii(2-3)此即Verlet算法的速度形式(A3),其特點(diǎn)是:(1)啟動(dòng)此算法，必須已知兩點(diǎn)r0,v0 ；i i(2) r n+1 和Vn+1 同步計(jì)算。(Fn+1的計(jì)算需要知道 rn+1 )此算法比A2算法更穩(wěn)定。3.趨恒階段的能量調(diào)整由于系統(tǒng)初始狀態(tài)的能量離我們要模擬系統(tǒng)的能量(或溫度

29、)有一定差 異，這就需要在系統(tǒng)趨于平衡的階段進(jìn)行調(diào)整。最常用的方法之一就是進(jìn) 行速度標(biāo)度：Vn+1 J & Vn+1標(biāo)度因子:p =，0N _ 3)kT /Z mv2i/2refiT為設(shè)定的系統(tǒng)的溫度值。ref能量設(shè)定值：E =(3N-3)kTref 2ref能量參考值： = lz mv22，E皿=P 2 EE = m(pv )2ref 2i要注意的問題：一般不需要每一步都進(jìn)行標(biāo)度，可每隔若干步(比如50步)調(diào)整一次；當(dāng)能量達(dá)到所設(shè)定的值后，停止標(biāo)度，在以下的模擬時(shí)間內(nèi)能量保持不變。4.實(shí)例氯原子系統(tǒng)，N=256LJ 勢(shì)： TOC o 1-5 h z u7 n作用力:Fx=8u(r.)8xi4

30、8-X(b)14(b)取時(shí)間單位:1/2取長(zhǎng)度單位:取質(zhì)量單位:Fx(rij)=483 x 10-12 sm=(8 )10-26 kg（氯原子質(zhì)量）(b)14(b)48814(bmb 2使用上述單位，可得到約化單位下的作用力表達(dá)式:*i-x*) jr *k ij標(biāo)度因子:p = (3N - 3)燈局 / E mv21/2=3(N - 1)kTref/mb2 m| Ev*21/2=(N - 1) kTr / 16Ev*21/2 8i=(N-1)T*/16E v*21/2i約化溫度:kAg AT* = kTf T /gg = 119.8Kk取時(shí)間步長(zhǎng)h = 2 10-14S20fs,相當(dāng)于約化時(shí)間

31、約化溫度：T* =303K, T* =約化數(shù)密度：* = N/L3, N = 256L =是錯(cuò)的，* = 0.636 L = ,* = 0.83134當(dāng) N = 64 時(shí)， * = 0.83134 L =,所以，取.=因?yàn)橐髍c L/2 =。(習(xí)題)計(jì)算源程序見p129,程序PL1正則系綜分子動(dòng)力學(xué)方法正則系綜：NVT恒定，且P = 0如何保持溫度T恒定1.速度標(biāo)度方法Vn+1 J & Vn+1p =(3N - 3)kT / Z mv2i/2、efi（1）規(guī)定初始位置見和初始速度v o ；（2）計(jì)算 J rn+1 = rn + *n + 2F n /2m ；v n+1 = v n + h（F

32、 n+1 + F n ） / 2m（3）計(jì)算 p =（3N - 3）kT /z mv2i/2 ；refi對(duì)速度進(jìn)行標(biāo)度： Vn+1 P Vn+1， 返回（2）。說明：（1）在前面的微正則系綜的模擬中也對(duì)速度進(jìn)行了標(biāo)度，但其目的 是使總能量達(dá)到所需要得值，而且不必每步進(jìn)行速度標(biāo)度，可 隔若干步標(biāo)度一次，一旦總能量達(dá)到所設(shè)定的值，即可停止標(biāo) 度；而在正則系綜的模擬中，每步都必須進(jìn)行標(biāo)度，以始終保 持溫度恒定。（2）通過一個(gè)廣義位勢(shì)引進(jìn)能量漲落，連同細(xì)致耦合的一種特殊選 擇，會(huì)導(dǎo)致速度標(biāo)度機(jī)制。（見書中的證明，公式有誤）（3）由于控制溫度的反饋環(huán)內(nèi)的時(shí)間延遲，溫度將在一定程度上漲 落。2.隨機(jī)方法（

33、Anderson熱?。w系與一強(qiáng)加了溫度的熱浴相耦合，與熱浴的耦合由偶爾作用于隨機(jī)選擇的粒子上的脈沖力表示。在開始模擬前首先要選擇系統(tǒng)與熱浴的耦合強(qiáng)度，這一耦合強(qiáng)度由 隨機(jī)碰撞的頻率 決定。模擬步驟:（1）規(guī)定初始位置見和初始速度W ；（2）計(jì)算rn+i = rn + hvn + h2Fn /2m .vn+1 = vn + h（ F n+1 + Fn） / 2m（3）粒子i在時(shí)間間隔h內(nèi)發(fā)生碰撞的幾率為h ，產(chǎn)生隨機(jī) 數(shù)，如果 h ,則粒子i經(jīng)歷一次隨即碰撞，否則，返回；（4）粒子i從溫度為T的麥克斯韋一玻爾茲曼分布中獲得一個(gè) 新的速度，返回（2）。等溫等壓系綜分子動(dòng)力學(xué)方法1,等壓等焓系綜N

34、 P H恒定，P = 0 （P：壓強(qiáng)，P：總動(dòng)量）H = E + PV（PE為外壓強(qiáng)，平衡時(shí)，Pe=P）AH = AE + pAV若系統(tǒng)與外界絕熱，則E = -罕,即曲二（等焓過程）要保持壓強(qiáng)不變，體積則必須變化，為此，將體積V作為一個(gè)新的動(dòng)力學(xué)變量，構(gòu)造一個(gè)新的拉氏函數(shù): .L(r, r,V,V) = T - U + T - U1 -=_ mr 2 - U (r) + MV 2 - PVi2 ei引入標(biāo)度坐標(biāo):p = r /V1/3 = r /L 或 r = Lpr（=u = p）= Lp + Lp夜Lp（ L 0，體積變化緩慢）. bl， v - mV (共軛動(dòng)量)p L mp =如iL(p,