計(jì)算物理第六章

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

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

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

4、則運(yùn)動(dòng)的Lagrangian方程為d ¶ L¶ L= 0t d單原子的運(yùn)動(dòng)方程p=r&iim=if i=-Ñmi &r&iVr i數(shù)值求解時(shí)的離散化方法對(duì)體系的運(yùn)動(dòng)方程組采用計(jì)算機(jī)進(jìn)行數(shù)值求解時(shí)，需要將運(yùn)動(dòng)方程離散化為有限差分方-程。常用的求解方法有歐、法等。數(shù)值計(jì)算的誤差階數(shù)顯然取決于所采用的數(shù)值求解方法的近似階數(shù)。原則上， 只要計(jì)算機(jī)計(jì)算速度足夠大，內(nèi)存足夠多，我 們可以使計(jì)算誤差足夠小。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é)動(dòng)力學(xué)模擬的基本步驟動(dòng)力學(xué)模擬的實(shí)際步驟可以劃分為四步：設(shè)定模擬所采用的模型；給定初始條件；趨于平衡的計(jì)算過程；宏 觀物理量的計(jì)算。1模擬模型的設(shè)定硬球勢(shì)ì +¥ ,0 ,如果r < s ,U (r ) =í如果r ³î.Lennard-Jones型勢(shì)位勢(shì)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根據(jù)經(jīng)典物理學(xué)的規(guī)律我們就可以知道在系綜模擬中的守恒量。微正則系綜的模擬中能量、動(dòng)量和角動(dòng)量均為守恒 量。在此系綜中他們分別表示為：vm(r ) + V (r )ùé1P = å piM = å r ´ på

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

9、應(yīng)當(dāng)選得很大。有限r(nóng)r - rr+ nrLr= min粒子間相互作用的最小像力約定,ijij動(dòng)力學(xué)模擬的最小像力約定示意圖采用最小像力約定后，元胞內(nèi)第 i 個(gè)粒子與周圍粒子的相互作用勢(shì)和相互作用力為Ui ( )=ri ¹ j( )åRu rijj =1, Nr (r)=i ¹ jå F ( )FiRrijrij.j =1, NR = rr , rr ,.rr,表示元胞內(nèi)所有粒子的坐標(biāo)。 rij是沿12Nrj - ri方向的矢量采用最小像力約定會(huì)使得在截?cái)嗵幜Ｗ拥氖芰τ幸粋€(gè)d -函數(shù)的奇異性，這會(huì)給模擬計(jì)算帶來誤差。為減小這種誤差，我們總可以將相互作用勢(shì)能

10、移到V (r ) - V (rc )，以保證在截?cái)嗵幭嗷プ饔媒咏恪?.給定初始條件給定粒子的初始位置和速度的數(shù)值：（1） 令初始位置在差分網(wǎng)格格子上，初始速度從玻爾茲曼分布隨機(jī)抽樣得到。（2） 令初始位置隨機(jī)地偏離差分網(wǎng)格格子，初始速度為零。（3） 令初始位置隨機(jī)地偏離差分網(wǎng)格格子，初始速度從玻爾茲曼分 布隨機(jī)抽樣得到。3趨于平衡使系統(tǒng)達(dá)到平衡，模擬中需要一趨衡過程。在這個(gè)過程中，我們?cè)黾踊驈南到y(tǒng)中移出能量，直到系統(tǒng)具有所要求的能量。4宏觀物理量的計(jì)算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)動(dòng)力學(xué)模擬1.微正則系綜的動(dòng)力學(xué)模擬粒子數(shù)恒定、體積恒定、能量恒定、整個(gè)系統(tǒng)的總動(dòng) 量恒等于零。(N,V,E)動(dòng)力學(xué)模擬步驟如下（Verlet算法）：（1）給定初始空間位置。（2）計(jì)算在第n

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

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