計算物理第六章

1、第六章動力學方法Molecular Dynamics Simulations第一節(jié)引言計算機模擬分類:(1)隨機模擬方法。優(yōu)點:隨機模擬方法計算的程序簡單，占內存少，但是該方法難于處理非平衡態(tài)的問題。該方法可以被用到 沒有任何內稟動力學模型體系的模擬上。(2)動力學方法（MolecularDynamics或簡稱MD）。可以處理非平衡態(tài)問題。但是使用該方法的程序較復 雜，計算量大，占內存也多。動力學方法利用相空間中的軌跡。素。古典力學來計算許多在動力學方法中不存在任何隨機因動力學簡介v動力學方法（MolecularDynamics或簡稱MD）是、水平上求解多體問題的重要的計算機模擬方法。v動力學

2、方法可以用于模擬與粒子運動路徑相關的基本過程。v在動力學中，粒子的運動行為是利用古典力學來計算許多原子、在相空間中的軌跡。系統(tǒng)的動力學機制決定運動方程的形式在動力學方法處理過程中，方程組的建立是通過對物理體系的微觀數學描述給出的。每個都各自服從經典的力學。這種方法可以處理與時間有關的過程，因而可以處理非平衡態(tài)問題。但是使用該方法的程序較復雜，計算量大，占內存也多。適用范圍廣泛原則上，動力學方法所適用的微觀物理體系并無什么限制。這個方法適用的體系既可以是少體系統(tǒng)，也可以是多體系統(tǒng)；既可以是點粒子體系，也可以是 具有內部結構的體系；處理的微觀客體既可以是也可以是其它的微觀粒子。，實際使用的限制實際

3、上，動力學模擬方法和隨機模擬方法一樣都著兩個基本限制：(1)有限觀測時間的限制；(2)有限系統(tǒng)大小的限制。通常人們感的是體系在熱力學極限下（即粒子數目趨于無窮時）的性質。但是計算機模擬的體系大小要比熱力學極限小得多，因此可能會出現有限 效應。為了減小有限效應，人們往往引入周期性、全反射、漫反射等邊界條件。當然邊界條件的引入顯然 會影響體系的某些性質。計算機平均值之計算1rrt()( )()òA = lim, p t=dtA r tAt ®¥ t0溫度之計算1 mv2= 1 k TB22第二節(jié)動力學基礎知識粒子運動方程及其數值求解：1.定義Lagrangian函數為

4、則運動的Lagrangian方程為d ¶ L¶ L= 0t d單原子的運動方程p=r&iim=if i=-Ñmi &r&iVr i數值求解時的離散化方法對體系的運動方程組采用計算機進行數值求解時，需要將運動方程離散化為有限差分方-程。常用的求解方法有歐、法等。數值計算的誤差階數顯然取決于所采用的數值求解方法的近似階數。原則上， 只要計算機計算速度足夠大，內存足夠多，我 們可以使計算誤差足夠小。Verlgorithmr(t+t) = r(t) + v(t)t + (1/2)a(t)t2(1)r(t-t)= r(t) v(t)t+ (1/2)a

5、(t)t2(2)將上面兩式相加，得到：r(t+t) = 2r(t) r(t- t) + a(t)t2(3)v(t+t) = v(t) + a(t)t + (1/2)b(t)t2 a(t+t) = a(t) + b(t)(4)將 (5) 式的b(t) 代入 (4) ，得到：v(t+t) = v(t) + (1/2)a(t) + a(t+t) t(6)其他求解算法：Leap-frog algorithmr(t+t)= r(t) + v(t+(1/2)t) tv(t+(1/2)t) = v(t-(1/2)t) + a(t) tBeemans algorithmr(t+t) = r(t) + v(t)

6、t + (2/3)a(t)t2 (1/6)a(t-t)t2v(t+t) = v(t) + v(t)t + (1/3)a(t)t +(5/6)a(t)t(1/6)a(tt)t第三節(jié)動力學模擬的基本步驟動力學模擬的實際步驟可以劃分為四步：設定模擬所采用的模型；給定初始條件；趨于平衡的計算過程；宏 觀物理量的計算。1模擬模型的設定硬球勢ì +¥ ,0 ,如果r < s ,U (r ) =í如果r ³î.Lennard-Jones型勢位勢V(r)43力F(r)2éæ s ö12æ s ö6 

7、49;1U (r ) = 4e êç- ç÷÷úúû0êëè r øè r ø-1-2排斥力吸引力 -30.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6v根據經典物理學的規(guī)律我們就可以知道在系綜模擬中的守恒量。微正則系綜的模擬中能量、動量和角動量均為守恒 量。在此系綜中他們分別表示為：vm(r ) + V (r )ùé1P = å piM = å r ´ på

8、2&E =êúiiiië2ûiii元胞周期性邊界條件,該邊界條件的具體實現是這樣操作的：當有一個粒子穿過基本動力學元胞的六方體表面時，就讓這個粒子以相同的速度穿過此表面對面的表面重新進入動力學元胞內。rc不同動力學元胞盒子內粒子間的相互作用：對于不同動力學元胞盒子內粒子間的相互作用，如果相互作用是短程力，我們可以在長度rc處截斷。這里V(rc)必須要足夠小，以使截斷顯著地影響模擬結果。型力學元胞尺度L通常選得比rc大很多。我們往往選擇元胞尺度滿足不等式條件L/2> rc ，使得距離大于L/2的粒子的相互作用可以忽略，以避免效應。通常L的數值

9、應當選得很大。有限rr - rr+ nrLr= min粒子間相互作用的最小像力約定,ijij動力學模擬的最小像力約定示意圖采用最小像力約定后，元胞內第 i 個粒子與周圍粒子的相互作用勢和相互作用力為Ui ( )=ri ¹ j( )åRu rijj =1, Nr (r)=i ¹ jå F ( )FiRrijrij.j =1, NR = rr , rr ,.rr,表示元胞內所有粒子的坐標。 rij是沿12Nrj - ri方向的矢量采用最小像力約定會使得在截斷處粒子的受力有一個d -函數的奇異性，這會給模擬計算帶來誤差。為減小這種誤差，我們總可以將相互作用勢能

10、移到V (r ) - V (rc )，以保證在截斷處相互作用接近零。2.給定初始條件給定粒子的初始位置和速度的數值：（1） 令初始位置在差分網格格子上，初始速度從玻爾茲曼分布隨機抽樣得到。（2） 令初始位置隨機地偏離差分網格格子，初始速度為零。（3） 令初始位置隨機地偏離差分網格格子，初始速度從玻爾茲曼分 布隨機抽樣得到。3趨于平衡使系統(tǒng)達到平衡，模擬中需要一趨衡過程。在這個過程中，我們增加或從系統(tǒng)中移出能量，直到系統(tǒng)具有所要求的能量。4宏觀物理量的計算A(rr ( N ) ( ), pr ( N ) ( ) 1(t ¢ - t0 )¢tòA = limdt

11、62;®¥t0dtE ( p (N ) (t )t¢r1(t ¢ - t0 )òtE= limkkt¢®¥0( p 2 )( m )1nNååiE=>kn - n2mm >n0 i=10u()n=1åå( m )Urijn - nm >n0 i< j0第四節(jié) 平衡態(tài)動力學模擬1.微正則系綜的動力學模擬粒子數恒定、體積恒定、能量恒定、整個系統(tǒng)的總動 量恒等于零。(N,V,E)動力學模擬步驟如下（Verlet算法）：（1）給定初始空間位置。（2）計算在第n

12、步時粒子所受的力。（3）計算在第n+1步時所有粒子所處的空間位置。r (n)2r (n+1) = 2r (n) - r (n-1)+ Fih/ miii（4）計算第n步的速度。=( n +1)- r( n -1)v ( n)r/ 2hiii（5）返回到步驟（2），開始的模擬計算。Verlet算法的速度形式：r(1)r（1）（2）給定初始空間位置。ir(1)v給定初始速度。i（3）利用公式：rr ( n+1) = r ( n) + hv ( n)+ Fih / 2m( n)2iii( n +1)計算在第n+1步時所有粒子所處的空間位置r。i（4）計算在第n+1步時所有粒子的速度：+ h(+ Fi

13、)/ 2mrrr( n +( n += v1)( n )1)( n )vFiii（5）返回到步驟（3），開始第n+2步的模擬計算。速度標度因子:ù1 / 2éb = êT (N - 1) ú*êúå i216vêúëûi2.正則系綜的動力學模擬ù1/ 2é速度標度因子:b = ê(3N - 4)kT úêúåi2mvêúëiû正則系綜動力學的模擬具體步驟:(Verlet算法的速度形式)（1）給定初始空間位置，（2）給定初始速度，（3）利用公式：rrr(n+= r+(n) + F (n) h 21)( n)rhv/ 2miiii計算在第n+1步時所有粒子所處的空間位置，（4） 計在時所有粒子的速度：+ h()/ 2mrr( n +( n += v+1)( n )1)( n )vFFiiii動能和速度標度因子：ù1 / 2&#