物理問題的計算機模擬方法(1)-分子動力學(xué)

物理問題的計算機模擬方法(1)—分子動力學(xué)編制僅供參考審核批準(zhǔn)生效日期地址：電話：傳真:郵編:碩士研究生課程《物理問題的計算機模擬方法》講義適用專業(yè):凝聚態(tài)物理、材料物理與化學(xué)、理論物理、光學(xué)工程學(xué)時：30—40學(xué)時參考教材：1．[德]D.W.Heermann著，秦克誠譯，理論物理中的計算機模擬方法，北京大學(xué)出版社，1996。2．[荷]Frenkel&Smit著，汪文川等譯，分子模擬—從算法到應(yīng)用，化學(xué)工業(yè)出版社，2002。3．M.P.AllenandD.J.Tildesley,ComputerSimulationofLiquids,ClarendonPress,Oxford,1989.4.A.R.Leach,MolecularModelling:PrinciplesandApplications,AddisonWesleyLongman,England,1996.5.[德]D.羅伯著，計算材料學(xué)，化學(xué)工業(yè)出版社，2002。6.[英]B.Chopard&MichelDroz著，物理系統(tǒng)的元胞自動機模擬，祝玉學(xué)，趙學(xué)龍譯，清華大學(xué)出版社，2003。

目錄計算機模擬方法概論1.1序言1.2熱力學(xué)系統(tǒng)物理量的統(tǒng)計平均1.3分子動力學(xué)方法模擬的基本思想1.4蒙特卡羅方法模擬的基本思想1.5元胞自動機模擬的基本思想1.5.1簡要的發(fā)展歷程1.5.2簡單元胞自動機：奇偶規(guī)則1.5.3元胞自動機的一般定義確定性模擬方法—分子動力學(xué)方法（MD）2.1分子動力學(xué)方法2.2微正則系綜分子動力學(xué)方法2.3正則系綜分子動力學(xué)方法2.4等溫等壓系綜分子動力學(xué)方法隨機性模擬方法—蒙特卡羅方法（MC）3.1預(yù)備知識3.2布朗動力學(xué)（BD）3.3蒙特卡羅方法3.4微正則系綜蒙特卡羅方法3.5正則系綜蒙特卡羅方法3.6等溫等壓系綜蒙特卡羅方法3.7巨正則系綜蒙特卡羅方法離散性模擬方法—原胞自動機（CA）4.1引言4.2元胞自動機模擬*4.3元胞自動機模擬的應(yīng)用

第一章計算機模擬方法概論§1.1序言什么是計算機模擬？SimulationModelling2．為什么要進行計算機模擬？3．常用的計算機模擬方法確定性模擬方法：MD模擬MolecularDynamics隨機性模擬方法：MC模擬MonteCarlo離散性模擬方法：CA模擬CellularAutomata§1.2熱力學(xué)系統(tǒng)物理量的統(tǒng)計平均描述系統(tǒng)的坐標(biāo)（自由度）：x(t)={x1(t),x2(t),…xN(t)}系統(tǒng)的物理量：A(x(t))1．時間平均←分子動力學(xué)（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)的平均，則有實際模擬的系統(tǒng)大小也是有限的：有限的粒子數(shù)N或有限的系統(tǒng)限度L對統(tǒng)計平均結(jié)果有影響。§1.3分子動力學(xué)(MD)方法模擬的基本思想基本原理系統(tǒng)：N個粒子，體積V，粒子質(zhì)量為m描述一個粒子運動狀態(tài)的自由度：（ri,pi）（pi=mvi）相空間：6N維，相空間中的一點的坐標(biāo)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)：實際模擬時間總是有限的，模擬時間的長短可通過判斷時間的增加對平均值的影響來確定，當(dāng)繼續(xù)增加時間帶來的平均值得變化在允許的誤差范圍之內(nèi)時，即可認(rèn)為模擬足夠長了。計算步驟運動方程：即(1-6)或(1-7)數(shù)值求解：用差分近似表示微分（采用不同的差分格式，可得到不同的算法）。用顯示中心差分格式，將（7）式寫為(1-8)由（7）和（8）式可得：(1-9)第一步：由（9）式計算第i個粒子在t+Δt時刻的位置坐標(biāo)。要啟動計算，我們必須要知道最初兩點ri(0)和ri(Δt)第二步：對不同時刻t=Δt,2Δt,3Δt,……,LΔt(t0=0)計算物理量A(r1(lΔt),r2(lΔt),……,rN(lΔt))(l=1,2,……,N)第三步：計算物理量A的平均值L的大小由繼續(xù)增大L而<A>不變（或變化在誤差范圍內(nèi)）來確定。§1.4蒙特卡羅（MC）方法模擬的基本思想基本原理以正則系綜（T,V,N）為例正則分布：正則配分函數(shù)：系統(tǒng)能量：物理量：A(rN)=A(r1,r1,…,rN)系綜平均： （1-10）（位形積分）（1-11）用MC方法計算上述多維積分。計算步驟劃分原胞N個粒子—3N個空間自由度，3N維空間劃分成s個相等的原胞（s>>1）注意：由于積分中不含動量，所以我們只需要在位置空間積分，而不需要在相空間中積分。當(dāng)系統(tǒng)的代表點落入第i個原胞時，則認(rèn)為系統(tǒng)處在狀態(tài)i，因此，s為系統(tǒng)可能的微觀狀態(tài)數(shù)目。于是，積分（10）和（11）可近似表述為(1-12)(1-13)建立馬爾可夫(Maρkoв)過程（鏈）將s個狀態(tài)看作一組隨機事件馬爾可夫鏈：從狀態(tài)i→j狀態(tài)j(i→j)的概率pij，只與i和j有關(guān)。，i=1,2,…,s若i經(jīng)歷n步到達j，其概率表示為，存在極限概率(j=1,2,…,s)uj為系統(tǒng)處在狀態(tài)j的概率。于是，沿?zé)o限長的馬爾可夫鏈，物理量A的平均值可寫為(1-14)選取，則（14）式為A的正則系綜平均值。（3）抽樣方法采用怎樣的抽樣方法所構(gòu)成的馬爾可夫鏈能得到上述平均值？粒子位置坐標(biāo)：粒子編號：=1,2,…N坐標(biāo)的三個方向：=1,2,3系統(tǒng)狀態(tài)：i=1,2,…,s給定粒子位置坐標(biāo)的變化量(小于系統(tǒng)體積的限度)給定系統(tǒng)的初態(tài)i，隨機選定4個隨機數(shù)，其中三個（=1，2，3），且–11，一個表示粒子編號=1,2,…N，由此隨機確定粒子位置的變化：（確保）若，則運動到新的位置，即系統(tǒng)由狀態(tài)i過渡到狀態(tài)j；若，則再選一個隨機數(shù)4（041），若，則粒子保留在原位置，不發(fā)生ij的躍遷；若，則發(fā)生ij的躍遷。由此進行下去，則形成一個馬爾可夫鏈（或過程），此鏈的長度L(即粒子行走的步數(shù)，遠(yuǎn)大于s)，由所計算的物理量的平均值(1-15)不再隨鏈的加長而改變來確定。由此得到的平均值即正則平均值。一般來說，L與N,V,T有關(guān)，比如，N=32~108,L=3000~5000。歸納起來，計算系統(tǒng)物理量的正則系綜平均值的具體步驟如下：第一步：給定系統(tǒng)的初始狀態(tài)(粒子的初始位置)ri和每一步的改變量；第二步：選擇四個隨機數(shù)，其中一個代表粒子的編號i(1iN)；另外三個表示粒子空間坐標(biāo)的改變x,y,z(-,=1,2,3);第三步：計算粒子i的新位置第四步：計算粒子在新舊兩個位置系統(tǒng)的能量之差第五步：由的大小判斷粒子i是否從ri運動到:若，則ri；若，則再選一個隨機數(shù)R(0<R<1),如果，則ri；如果，則ri不變，返回第二步。第六步：計算第七步：重復(fù)上述各個步驟，直到完成L步為止，最后利用公式（15）計算A的平均值。粒子間相互作用勢模型的選取最簡單的兩種模型：（1）硬球模型(為硬球的直徑)（2）L—J勢§1.5元胞自動機(CA)模擬的基本思想元胞自動機：時間和空間都離散、物理參量只取有限數(shù)值集的物理系統(tǒng)的理想化模型cellularautomata或cellularautomaton—CA1.5.1簡要的發(fā)展歷程1．自繁殖系統(tǒng)20世紀(jì)40年代，VonNeumann,構(gòu)造能解決非常復(fù)雜問題的計算機，設(shè)想模仿人腦的行為——尋求與生物過程無關(guān)的情況下自繁殖機理的邏輯抽象。根據(jù)S.Ulam的建議，VonNeumann在由元胞構(gòu)成的完全離域的框架下處理這個問題，構(gòu)造了一個完全離散的動力學(xué)系統(tǒng)——元胞自動機。第一個自復(fù)制元胞自動機——由二維方形網(wǎng)絡(luò)組成，有數(shù)千個基本元胞構(gòu)成的自繁殖結(jié)構(gòu)。（1）一般機器只能構(gòu)造比自己簡單的客體，而采用自復(fù)制元胞機，可獲得一種能產(chǎn)生新的、具有同樣復(fù)雜性和功能的“機器”；（2）VonNeumann的元胞自動機規(guī)則具有所謂通用計算的性質(zhì)，這意味著，存在一種元胞自動機的初始構(gòu)型，該元胞自動機能產(chǎn)生任何計算機算法的解。通用計算的性質(zhì)指：用元胞自動機演化規(guī)則能夠模擬任何計算機流程（邏輯選擇器開關(guān)）。2．生命游戲機1970年，數(shù)學(xué)家JohnConway生命游戲機的概念，尋找能導(dǎo)致復(fù)雜行為的簡單規(guī)則。設(shè)想一個類似于棋盤的二維方形網(wǎng)格，每個元胞可能的狀態(tài)是活（狀態(tài)1）或死（狀態(tài)0），其更新規(guī)則是：有三個活元胞包圍的一個死元胞恢復(fù)為活元胞；由兩個以下或三個以上活元胞包圍的活元胞因孤立或擁擠而死亡。結(jié)果表明，生命游戲機有出乎意料的豐富行為，從原“湯”中顯示出來的復(fù)雜結(jié)構(gòu)，演變發(fā)展成為某些特殊的技藝，例如，可能形成所謂的滑翔機——緊鄰元胞的特殊排列，這些元胞具有沿直線彈道穿越空間運動的特性。生命游戲機也是具有計算通用性的元胞自動機。3．模擬物理系統(tǒng)（1）20世紀(jì)70年代，Hardy,Pomeau和Pazzis建立了所謂的HPP格子氣體模型，用以在質(zhì)量和動量守恒的情況下在方形網(wǎng)格上模擬粒子的碰撞行為。（2）1986年，F(xiàn)risch,Hasslacher和Pomeau提出了著名的FHP模型，這是在六邊形網(wǎng)格上模擬二維流體動力學(xué)的第一個嚴(yán)格模型——全離散計算機模型替代風(fēng)洞試驗。HPP和FHP通常稱之為格子氣自動機（LGA—LatticeGasAutomata）（3）Ising自旋動力學(xué)模型，20世紀(jì)80年代末，Vichniac提出Q2R規(guī)則。（4）格子Boltzmann方法與多粒子模型格子Boltzmann方法或模型（LBM）：網(wǎng)格上定義的物理模型，在這個網(wǎng)格上與每個格位相關(guān)聯(lián)的變量是平均粒子數(shù)，或具有一定速度粒子出現(xiàn)的概率。該模型可以用均化或因子分解方法由元胞自動機動力學(xué)推導(dǎo)出來，或者自定義，而與特定的實現(xiàn)無關(guān)。格子Boltzmann模型保持了元胞自動機方法的微觀水平解釋，但忽略了多體的相關(guān)函數(shù)，但這種方法已經(jīng)成為目前模擬物理系統(tǒng)中最有前途的方法之一。在嚴(yán)格的元胞自動機方法與較靈活的格子Boltzmann模型之間，有一種目前處于發(fā)展中的模型——多粒子模型。這種模型保留量化狀態(tài)的概念，但接受無限數(shù)值集，因此既保證了數(shù)值穩(wěn)定性（與LBM相反），又考慮了多體相關(guān)性。在模擬物理系統(tǒng)時，大量的可能狀態(tài)提供更多的靈活性，并產(chǎn)生小的統(tǒng)計噪聲。但多粒子動力學(xué)更難設(shè)計，且在數(shù)值計算上比格子Boltzmann方法更慢。1.5.2簡單元胞自動機：奇偶規(guī)則簡單元胞自動機的演化規(guī)則（20世紀(jì)70年代，EdwardFredkin提出，定義在二維方形網(wǎng)上）：網(wǎng)格的每一個格位是一個元胞，以其位置r=(i,j)來標(biāo)記，其中i和j為行和列的標(biāo)號。函數(shù)描述每個元胞在時間t的狀態(tài)，其值為0或1。從時間t=0的初始條件及網(wǎng)格上給定的構(gòu)形值開始，t=1時狀態(tài)按下列步驟求得：（1）對于每個格位r，都計算出其位于東、南、西、北4個最近鄰格位的值之和。應(yīng)使系統(tǒng)在i和j兩個方向循環(huán)（如同在環(huán)面上），從而確定出所有格位的計算值。（2）如果這個和值為偶數(shù)，則新狀態(tài)為0（白色），否則為1（黑色）。重復(fù)上述步驟，得到t=2,3,4,…的狀態(tài)。這個元胞自動機的奇偶規(guī)則可表示為：式中，符號代表“異或”邏輯運算，也即模2和：11=00=0，10=01=1。反復(fù)迭代這個規(guī)則時，可得到非常精美的幾何圖形。1.5.3元胞自動機的一般定義定義：規(guī)整的元胞網(wǎng)格覆蓋d維空間的一部分；歸屬于網(wǎng)格的每個格位r的一組布爾變量(r,t)={1(r,t),2(r,t),…,m(r,t)}給出每個元胞在時間t=0,1,2,…的局部狀態(tài)；演化規(guī)則R={R1,R2,…,Rm}按下列方式指定狀態(tài)(r,t)的時間演化過程：j(r,t+1)=Rj[(r,t),(r+1,t),(r+2,t),…,(r+q,t)]式中r+q指定從屬于元胞r的給定鄰居。鄰居：二維元胞自動機的兩種鄰居：（1）VonNeumann鄰居，有一個中心元胞（要演化的元胞）和四個位于其近鄰東西南北方位的元胞組成；（2）Moore鄰居，除了前面涉及的最近鄰元胞外，還包括次近鄰的4個元胞，共九個元胞。還有一種有用的鄰居稱為Margolus鄰居，將空間劃分成22元胞的鄰接單元塊，這個規(guī)則對位于單元塊內(nèi)的位置即左上、右上和左下、右下很敏感。三種鄰居如下圖所示。邊界條件：周期性邊界條件；固定邊界條件；絕熱邊界條件；映射邊界條件。備注：確定型元胞自動機：演化規(guī)則確定，給定的初始狀態(tài)將始終演化出同樣的式樣。概率型自動元胞機：演化規(guī)則包含一定的隨機性，給定的初始狀態(tài)可能演化出不同的式樣。

第二章確定性模擬方法—分子動力學(xué)方法（MD）§2.1分子動力學(xué)方法用分子動力學(xué)方法模擬氣體、液體或固體系統(tǒng)1．分子動力學(xué)的描述形式：哈密頓描述拉格朗日描述牛頓運動方程描述2．劃分MD元胞選取適當(dāng)大小的V=L3，太大耗費計算時間，太小不能準(zhǔn)確反映系統(tǒng)的性質(zhì)。目的：維持系統(tǒng)恒定的密度。對平衡態(tài)系統(tǒng)，液體和氣體原胞的形狀無關(guān)緊要，但對晶體則不然。周期性邊界條件目的：減少表面效應(yīng)相互作用勢、最小影像約定和切斷距離n=(n1,n2,n3)元胞內(nèi)的粒子間元胞內(nèi)的粒子與影像的相互作用粒子的相互作用最小影像約定(minimumimageconvention):最小影像元胞中的一個粒子只與元胞中的另外N-1個粒子中的每一個或其最近鄰影像發(fā)生相互作用。切斷距離（cutoffdistance）：rcL/2（切斷半徑）上述約定相當(dāng)于用rc來切斷位勢，即rij>rc的相互作用忽略不計，代價是忽略了背景。積分格式于是，運動方程可寫為(2-1)（1）遞推公式在第一章中求解運動方程時，我們是直接求解的關(guān)于粒子位置坐標(biāo)的二階微分方程，得到的遞推公式需要知道最初的兩點位置才能啟動計算，但在實際計算中，我們常常是給出最初的位置和速度，于是，我們可通過選取一定的差分格式，有運動方程(16)得到關(guān)于位置和速度的遞推公式。或已知遞推公式的具體形式取決于差分各式的選取。（2）時間步長選擇時間步長的原則：在保證計算精度的前提下，盡量節(jié)省計算時間。實例：Ar原子系統(tǒng)，L—J勢時間步長取為(=100ps=10fs)（3）約束條件保證能量、動量或角動量守恒。（4）減小計算誤差的技巧數(shù)值計算不可避免有誤差，與誤差有關(guān)的因素主要有：差分格式時間步長切斷半徑最小影像約定等等6.計算熱力學(xué)量若系統(tǒng)的NVE恒定，則（微正則系綜平均—microcannonicalensembleaverage）若系統(tǒng)的NVT恒定，則（正則系綜平均—cannonicalensembleaverage）在實際計算中，往往還需要要求系統(tǒng)的總動量P=0或恒定，因整個系統(tǒng)不受外界的作用。溫度與能量均分定理N—總粒子數(shù)，Nc—約束條件個數(shù),一般情況下，N>>NcP=0Nc=3位勢切斷誤差與修正g(r)—對關(guān)聯(lián)函數(shù)（paircorrelationfunction）—粒子數(shù)密度—當(dāng)原點位置上一個粒子時，在r附近dr內(nèi)發(fā)現(xiàn)有一個粒子的幾率。對氣體或液體，（各向同性）平均每個粒子的能量切斷誤差修正為：計算機模擬的組織初始化——給定粒子的初始位置、初始速度等；趨衡——從初始轉(zhuǎn)臺開始，按運動方程要求的規(guī)律，從非平衡態(tài)趨向平衡態(tài)；投產(chǎn)——計算物理量的統(tǒng)計平均值。§2.2微正則系綜分子動力學(xué)方法微正則系綜：NVE恒定，且P=0（分子動力學(xué)方法的自然選擇）1．Verlet算法運動方程的顯示中心差分格式由可得到如下遞推公式（2-2）此即Verlet算法（A2）,其特點是：要啟動此算法，必須已知兩點，所以又稱為二步法；和不能同步算出，第n+1步算出的第n步的速度。2．Verlet算法的速度形式由此可得如下遞推公式：（2-3）此即Verlet算法的速度形式（A3）,其特點是：啟動此算法，必須已知兩點；（2）和同步計算。（的計算需要知道）此算法比A2算法更穩(wěn)定。3．趨恒階段的能量調(diào)整由于系統(tǒng)初始狀態(tài)的能量離我們要模擬系統(tǒng)的能量（或溫度）有一定差異，這就需要在系統(tǒng)趨于平衡的階段進行調(diào)整。最常用的方法之一就是進行速度標(biāo)度：標(biāo)度因子：（2-4）Tref為設(shè)定的系統(tǒng)的溫度值。能量設(shè)定值：能量參考值：要注意的問題：一般不需要每一步都進行標(biāo)度，可每隔若干步（比如50步）調(diào)整一次；當(dāng)能量達到所設(shè)定的值后，停止標(biāo)度，在以下的模擬時間內(nèi)能量保持不變。實例氬原子系統(tǒng)，N=256

L作用力：取時間單位：取長度單位：=0.3405nm取質(zhì)量單位：m=6.6338210-26kg使用上述單位，可得到約化單位下的作用力表達式：標(biāo)度因子：約化溫度：取時間步長h=210-14s=20fs,相當(dāng)于約化時間0.064約化溫度：T*=2.53303K,T*=0.72286.5K約化數(shù)密度：*=N/L3,N=256*=0.636L=7.83,*=0.83134當(dāng)N=64時，*=0.83134L=4.25,所以，取rc=2.5是錯的，因為要求計算源程序見p129,程序PL1§2.3正則系綜分子動力學(xué)方法正則系綜：NVT恒定，且P=0如何保持溫度T恒定？1．速度標(biāo)度方法算法A4:規(guī)定初始位置和初始速度；計算；計算；對速度進行標(biāo)度：，返回（2）。說明：在前面的微正則系綜的模擬中也對速度進行了標(biāo)度，但其目的是使總能量達到所需要得值，而且不必每步進行速度標(biāo)度，可隔若干步標(biāo)度一次，一旦總能量達到所設(shè)定的值，即可停止標(biāo)度；而在正則系綜的模擬中，每步都必須進行標(biāo)度，以始終保持溫度恒定。通過一個廣義位勢引進能量漲落，連同細(xì)致耦合的一種特殊選擇，會導(dǎo)致速度標(biāo)度機制。（見書中的證明，公式3.59有誤）由于控制溫度的反饋環(huán)內(nèi)的時間延遲，溫度將在一定程度上漲落。2．隨機方法（Anderson熱?。w系與一強加了溫度的熱浴相耦合，與熱浴的耦合由偶爾作用于隨