1、碩士研究生課程物理問題的計算機模擬方法講義適用專業(yè): 凝聚態(tài)物理、材料物理與化學、理論物理、光學工程學 時：3040學時參考教材： 1德D.W.Heermann 著，秦克誠譯，理論物理中的計算機模擬方法，北京大學出版社，1996。 2荷 Frenkel & Smit 著，汪文川 等譯，分子模擬從算法到應用，化學工業(yè)出版社，2002。 3M.P.Allen and D.J.Tildesley, Computer Simulation of Liquids, Clarendon Press, Oxford, 1989. 4. A.R.Leach, Molecular Modelling:

2、Principles and Applications, Addison Wesley Longman, England, 1996. 5. 德 D.羅伯 著，計算材料學，化學工業(yè)出版社，2002。 6. 英 B. Chopard & Michel Droz 著，物理系統(tǒng)的元胞自動機模擬，祝玉學，趙學龍 譯，清華大學出版社，2003。目 錄第一章 計算機模擬方法概論1.1 序言1.2 熱力學系統(tǒng)物理量的統(tǒng)計平均1.3 分子動力學方法模擬的基本思想1.4 蒙特卡羅方法模擬的基本思想1.5 元胞自動機模擬的基本思想 簡要的發(fā)展歷程 簡單元胞自動機：奇偶規(guī)則 元胞自動機的一般定義第二章 確定

3、性模擬方法分子動力學方法（MD）2.1 分子動力學方法2.2 微正則系綜分子動力學方法2.3 正則系綜分子動力學方法2.4 等溫等壓系綜分子動力學方法第三章 隨機性模擬方法蒙特卡羅方法（MC）3.1 預備知識3.2 布朗動力學（BD）3.3 蒙特卡羅方法3.4 微正則系綜蒙特卡羅方法3.5 正則系綜蒙特卡羅方法3.6 等溫等壓系綜蒙特卡羅方法3.7 巨正則系綜蒙特卡羅方法第四章 離散性模擬方法原胞自動機（CA）4.1 引言4.2 元胞自動機模擬*4.3元胞自動機模擬的應用第一章 計算機模擬方法概論§ 1.1 序言1 什么是計算機模擬？ Simulation Modelling2為什么

4、要進行計算機模擬？3常用的計算機模擬方法確定性模擬方法：MD模擬 Molecular Dynamics隨機性模擬方法：MC模擬 Monte Carlo離散性模擬方法：CA模擬 Cellular Automata§ 1.2 熱力學系統(tǒng)物理量的統(tǒng)計平均描述系統(tǒng)的坐標（自由度）：x(t)=x1(t),x2(t),xN(t)系統(tǒng)的物理量：A(x(t)1時間平均 分子動力學（MD）模擬 (1-1)2系綜平均 蒙特卡羅(MC)模擬 (1-2) 分布函數(shù)（幾率密度函數(shù)） (1-3) 配分函數(shù) (1-4) 相空間 H(x)系統(tǒng)的哈密頓函數(shù)對于處于平衡態(tài)的系統(tǒng)，可以證明： 對于實際的有限時間內(nèi)的平均，

5、則有 實際模擬的系統(tǒng)大小也是有限的：有限的粒子數(shù)N或有限的系統(tǒng)限度L對統(tǒng)計平均結果有影響。§ 1.3 分子動力學(MD)方法模擬的基本思想1 基本原理系統(tǒng)：N個粒子，體積V，粒子質(zhì)量為m描述一個粒子運動狀態(tài)的自由度：（ri, pi） （pi=mvi）相空間：6N維，相空間中的一點的坐標 XN=rN, (mvN) rN=(r1, r2, , rN), vN=(v1, v2, , vN)粒子間的相互作用勢：U(rN)=U(r1, r2, , rN)=決定系統(tǒng)相軌跡XN（t）的運動方程： (1-5)物理量A的宏觀值，由A(XN) 的時間平均獲得，即 （離散情況：）對于平衡態(tài)： 實際模擬時間

6、總是有限的，模擬時間的長短可通過判斷時間的增加對平均值的影響來確定，當繼續(xù)增加時間帶來的平均值得變化在允許的誤差范圍之內(nèi)時，即可認為模擬足夠長了。2 計算步驟運動方程： 即 (1-6) 或 (1-7)數(shù)值求解：用差分近似表示微分 （采用不同的差分格式，可得到不同的算法）。用顯示中心差分格式，將（7）式寫為 (1-8)由（7）和（8）式可得： (1-9)第一步：由（9）式計算第i個粒子在t+t時刻的位置坐標。 要啟動計算，我們必須要知道最初兩點ri(0)和ri(t)第二步：對不同時刻 t = t , 2t, 3t, , Lt (t0 = 0) 計算物理量A(r1(lt), r2(lt), , r

7、N(lt) (l = 1, 2, , N)第三步：計算物理量A的平均值 L的大小由繼續(xù)增大L 而<A>不變（或變化在誤差范圍內(nèi)）來確定。§ 1.4 蒙特卡羅（MC）方法模擬的基本思想1 基本原理以正則系綜（T, V, N）為例正則分布：正則配分函數(shù)：系統(tǒng)能量：物理量：A(rN) = A(r1, r1, , rN )系綜平均： （1-10） （位形積分） （1-11）用MC方法計算上述多維積分。2 計算步驟（1） 劃分原胞 N個粒子3N個空間自由度，3N維空間劃分成s個相等的原胞（s>>1）注意：由于積分中不含動量，所以我們只需要在位置空間積分，而不需要在相空間

8、中積分。當系統(tǒng)的代表點落入第i個原胞時，則認為系統(tǒng)處在狀態(tài)i，因此，s為系統(tǒng)可能的微觀狀態(tài)數(shù)目。于是，積分（10）和（11）可近似表述為 (1-12) (1-13)（2） 建立馬爾可夫(Mako)過程（鏈）將s個狀態(tài)看作一組隨機事件馬爾可夫鏈：從狀態(tài)ij狀態(tài)j(eiej)的概率pij ，只與ei和ej 有關。 ，i=1, 2, , s若ei 經(jīng)歷n步到達ej ，其概率表示為，存在極限概率 ( j=1, 2, , s)uj為系統(tǒng)處在狀態(tài)j的概率。于是，沿無限長的馬爾可夫鏈，物理量A的平均值可寫為 (1-14) 選取 ，則（14）式為A的正則系綜平均值。（3）抽樣方法 采用怎樣的抽樣方法所構成的馬

9、爾可夫鏈能得到上述平均值？ 粒子位置坐標：粒子編號：g =1, 2, N坐標的三個方向：a = 1, 2, 3 系統(tǒng)狀態(tài)：i = 1, 2, , s給定粒子位置坐標的變化量d (小于系統(tǒng)體積的限度)給定系統(tǒng)的初態(tài)i，隨機選定4個隨機數(shù)，其中三個xa （a =1， 2，3），且 1 £ xa £1，一個表示粒子編號 g =1, 2, N，由此隨機確定粒子位置的變化： （ 確保 ）若，則 g 運動到新的位置，即系統(tǒng)由狀態(tài)i過渡到狀態(tài)j；若，則再選一個隨機數(shù)x 4（0 £ x 4 £1），若，則粒子保留在原位置，不發(fā)生i®j的躍遷；若，則發(fā)生i &#

10、174; j的躍遷。由此進行下去，則形成一個馬爾可夫鏈（或過程），此鏈的長度L(即粒子行走的步數(shù)，遠大于s)，由所計算的物理量的平均值 (1-15)不再隨鏈的加長而改變來確定。由此得到的平均值即正則平均值。一般來說，L與N, V, T 有關，比如，N=32108, L=30005000。歸納起來，計算系統(tǒng)物理量的正則系綜平均值的具體步驟如下：第一步：給定系統(tǒng)的初始狀態(tài)(粒子的初始位置)ri和每一步的改變量d；第二步：選擇四個隨機數(shù)，其中一個代表粒子的編號i ( 1£ i £ N )；另外三個表示粒子空間坐標的改變dx, dy, dz ( -d £ da £

11、;d , a=1, 2, 3);第三步：計算粒子i的新位置 第四步：計算粒子在新舊兩個位置系統(tǒng)的能量之差 第五步：由的大小判斷粒子i是否從ri運動到: 若，則ri®； 若，則再選一個隨機數(shù)R(0 < R < 1), 如果 ，則ri®； 如果 ，則ri 不變，返回第二步。第六步：計算第七步：重復上述各個步驟，直到完成L步為止，最后利用公式（15）計算A的平均值。3 粒子間相互作用勢模型的選取最簡單的兩種模型： （1）硬球模型 (s 為硬球的直徑) （2）LJ 勢§ 1.5 元胞自動機(CA)模擬的基本思想元胞自動機：時間和空間都離散、物理參量只取有限數(shù)值

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

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

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

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

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

17、= 1 時狀態(tài)按下列步驟求得：（1）對于每個格位r，都計算出其位于東、南、西、北4個最近鄰格位的值之和。應使系統(tǒng)在i和 j 兩個方向循環(huán)（如同在環(huán)面上），從而確定出所有格位的計算值。（2）如果這個和值為偶數(shù)，則新狀態(tài)為0（白色），否則為1（黑色）。重復上述步驟，得到t = 2, 3, 4, 的狀態(tài)。這個元胞自動機的奇偶規(guī)則可表示為：式中，符號Å代表“異或”邏輯運算，也即模2和：1Å1=0Å0=0，1Å0=0Å1=1。 反復迭代這個規(guī)則時，可得到非常精美的幾何圖形。 元胞自動機的一般定義定義：（1） 規(guī)整的元胞網(wǎng)格覆蓋d 維空間的一部分；（2）

18、歸屬于網(wǎng)格的每個格位r 的一組布爾變量F(r, t) = F1(r, t), F2(r, t), , Fm(r, t)給出每個元胞在時間t = 0, 1, 2, 的局部狀態(tài)；（3） 演化規(guī)則 R = R1, R2, , Rm按下列方式指定狀態(tài)F(r, t)的時間演化過程：Fj(r, t+1) = RjF(r, t), F(r + d1, t), F(r + d2, t), F(r + dq, t)式中r + dq 指定從屬于元胞r的給定鄰居。鄰居：二維元胞自動機的兩種鄰居：（1）Von Neumann鄰居，有一個中心元胞（要演化的元胞）和四個位于其近鄰東西南北方位的元胞組成；（2）Moore鄰

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

20、氣體、液體或固體系統(tǒng)1分子動力學的描述形式： 哈密頓描述 拉格朗日描述 牛頓運動方程描述2劃分MD元胞選取適當大小的 V=L3，太大耗費計算時間，太小不能準確反映系統(tǒng)的性質(zhì)。目的：維持系統(tǒng)恒定的密度。對平衡態(tài)系統(tǒng)，液體和氣體原胞的形狀無關緊要，但對晶體則不然。3 周期性邊界條件目的：減少表面效應4 相互作用勢、最小影像約定和切斷距離 n = (n1, n2, n3 ) ß ß 元胞內(nèi)的粒子間 元胞內(nèi)的粒子與影像的相互作用 粒子的相互作用最小影像約定(minimum image convention): 最小影像 元胞中的一個粒子只與元胞中的另外N-1個粒子中的每一個或其最近

21、鄰影像發(fā)生相互作用。切斷距離（cutoff distance）：rc £ L/2 （切斷半徑）上述約定相當于用rc 來切斷位勢，即rij > rc 的相互作用忽略不計，代價是忽略了背景。5 積分格式 于是，運動方程可寫為 (2-1) （1）遞推公式在第一章中求解運動方程時，我們是直接求解的關于粒子位置坐標的二階微分方程，得到的遞推公式需要知道最初的兩點位置才能啟動計算，但在實際計算中，我們常常是給出最初的位置和速度，于是，我們可通過選取一定的差分格式，有運動方程(16)得到關于位置和速度的遞推公式。 或 已知遞推公式的具體形式取決于差分各式的選取。（2）時間步長選擇時間步長的原

22、則：在保證計算精度的前提下，盡量節(jié)省計算時間。實例：Ar原子系統(tǒng)，LJ勢 時間步長取為 (=100ps=10fs)（3）約束條件保證能量、動量或角動量守恒。（4）減小計算誤差的技巧數(shù)值計算不可避免有誤差，與誤差有關的因素主要有： 差分格式 時間步長 切斷半徑 最小影像約定 等等6. 計算熱力學量 （1） 若系統(tǒng)的NVE恒定，則 （微正則系綜平均microcannonical ensemble average）（2） 若系統(tǒng)的NVT恒定，則 （正則系綜平均cannonical ensemble average）在實際計算中，往往還需要要求系統(tǒng)的總動量P=0或恒定，因整個系統(tǒng)不受外界的作用。（3）

23、 溫度與能量均分定理 N總粒子數(shù)，Nc約束條件個數(shù), 一般情況下，N>> Nc P=0 Þ Nc = 3（4） 位勢切斷誤差與修正g(r) 對關聯(lián)函數(shù)（pair correlation function） 粒子數(shù)密度當原點位置上一個粒子時，在r 附近dr內(nèi)發(fā)現(xiàn)有一個粒子的幾率。對氣體或液體， （各向同性）平均每個粒子的能量切斷誤差修正為：6 計算機模擬的組織（1） 初始化給定粒子的初始位置、初始速度等；（2） 趨衡從初始轉臺開始，按運動方程要求的規(guī)律，從非平衡態(tài)趨向平衡態(tài)；（3） 投產(chǎn)計算物理量的統(tǒng)計平均值。§ 2.2 微正則系綜分子動力學方法微正則系綜：NVE

24、恒定，且P = 0 （分子動力學方法的自然選擇）1Verlet 算法運動方程的顯示中心差分格式由可得到如下遞推公式 （2-2）此即Verlet 算法（A2）, 其特點是：（1） 要啟動此算法，必須已知兩點，所以又稱為二步法；（2） 和不能同步算出，第n+1步算出的第n步的速度。 2Verlet 算法的速度形式 由此可得如下遞推公式： （2-3）此即Verlet 算法的速度形式（A3）, 其特點是：（1） 啟動此算法，必須已知兩點；（2）和同步計算。（的計算需要知道）此算法比A2算法更穩(wěn)定。3趨恒階段的能量調(diào)整由于系統(tǒng)初始狀態(tài)的能量離我們要模擬系統(tǒng)的能量（或溫度）有一定差異，這就需要在系統(tǒng)趨于平

25、衡的階段進行調(diào)整。最常用的方法之一就是進行速度標度： 標度因子： （2-4）Tref 為設定的系統(tǒng)的溫度值。 能量設定值： 能量參考值： 要注意的問題：（1） 一般不需要每一步都進行標度，可每隔若干步（比如50步）調(diào)整一次；（2） 當能量達到所設定的值后，停止標度，在以下的模擬時間內(nèi)能量保持不變。4 實例氬原子系統(tǒng)， N=256LJ勢:作用力： 取時間單位：取長度單位：s = 0.3405nm取質(zhì)量單位：m=6.63382´10-26 kg(氬原子質(zhì)量)使用上述單位，可得到約化單位下的作用力表達式：標度因子： 約化溫度： 取時間步長h = 2´10-14s = 20fs,

26、相當于約化時間0.064約化溫度：T* = 2.53 ® 303K, T* = 0.722 ® 86.5K約化數(shù)密度：r* = N/L3, N = 256r* = 0.636 ® L = 7.83 , r* = 0.83134 ® L = 6.75當 N = 64 時，r* = 0.83134 ® L = 4.25, 所以，取rc = 2.5 是錯的，因為要求 rc > L/2 = 2.13。（習題3.5）計算源程序見p129, 程序PL1§ 2.3 正則系綜分子動力學方法正則系綜：NVT恒定，且P = 0 如何保持溫度T恒定？1速度標度方法 算法A4:(1) 規(guī)定初始位置和初始速度；(2) 計算；(3) 計算；(4) 對速度進行標度：，返回（2）。 說明：（1） 在前面的微正則系綜的模擬中也對速度進行了標度，但其目的是使總能量達到所需要得值，而且不必每步進行速度標度，可隔若干步標度一次，一旦總能量達到所設定的值，即可停止標度；而在正則系綜的模擬中，每步都必須進行標度，以始終保持溫度恒定。（2） 通過一個廣義位勢引進能量漲落，連同細致耦合的一種特殊選擇，會導致速度標度機制。（見書中的證明，公式3.59有誤）（3） 由于控制溫度的