3-分子動(dòng)力學(xué)方法-微正則系綜解析

三、分子動(dòng)力學(xué)方法

賴文生清華大學(xué)材料系技科樓2315室8/16/20241§3.1常規(guī)等容方法(NVE系綜)

分子動(dòng)力學(xué)模擬是一種用來計(jì)算一個(gè)經(jīng)典多體體系的平衡和傳遞性質(zhì)的方法。經(jīng)典這個(gè)詞意味著組成粒子的核心運(yùn)動(dòng)遵守經(jīng)典力學(xué)定律。真實(shí)實(shí)驗(yàn)：準(zhǔn)備試樣，將試樣同一測(cè)試儀器連結(jié)，在一定的時(shí)間間隔內(nèi)測(cè)量一下感興趣的性質(zhì)。如果測(cè)定服從統(tǒng)計(jì)噪音，那么測(cè)定的次數(shù)越多，結(jié)果就越準(zhǔn)確。分子動(dòng)力學(xué)模擬：準(zhǔn)備試樣：選擇一個(gè)由N個(gè)粒子組成的模型體系，解體系的牛頓運(yùn)動(dòng)方程直到體系的性質(zhì)不再隨時(shí)間改變（體系達(dá)到平衡）。平衡以后，就進(jìn)行實(shí)際的測(cè)定。模擬和實(shí)驗(yàn)容易犯的錯(cuò)誤很相似：如試樣未準(zhǔn)備好，測(cè)定時(shí)間太短，體系在實(shí)驗(yàn)過程中經(jīng)歷一不可逆的變化，或者并未測(cè)定相要測(cè)定的量。8/16/20242溫度在分子動(dòng)力學(xué)模擬中，測(cè)定一個(gè)可觀測(cè)的量，必須能將這一量表達(dá)為體系中粒子的位置與動(dòng)量的函數(shù)。例如，對(duì)溫度的方便定義是采用在所有自由度上的等配分能，以使它們以二次型進(jìn)入哈密爾頓函數(shù)。實(shí)際過程中，溫度是由測(cè)量體系的總動(dòng)能除以體系的自由度數(shù)Nf(對(duì)總動(dòng)能固定的體系Nf=3N-3，N是粒子數(shù))得到。由于體系總動(dòng)能漲落,體系的瞬時(shí)溫度也隨之改變溫度的相對(duì)漲落的數(shù)量級(jí)在左右。因此要達(dá)到溫度的一個(gè)準(zhǔn)確的估計(jì),應(yīng)對(duì)許多漲落進(jìn)行平均才行。

8/16/20243程序?qū)Ψ肿觿?dòng)力學(xué)模擬的最好介紹就是考察一個(gè)程序。程序按以下方式構(gòu)成。1）讀入指定運(yùn)算條件的參數(shù)（如初始溫度、粒子數(shù)、密度、時(shí)間等）。2）體系初始化（即選定初始位置和速度）。3）計(jì)算作用于所有粒子上的力4）解牛頓運(yùn)動(dòng)方程。這一步和上一步構(gòu)成了模擬的核心。重復(fù)這兩步直至我們計(jì)算體系的演化到指定的時(shí)間長(zhǎng)度。5）中心循環(huán)完成之后，計(jì)算并打印測(cè)定的平均值，模擬結(jié)束。8/16/20244虛擬算法Programmd

簡(jiǎn)單的MD程序callinit初始化t=0dowhile(t.lt.tmax)MD循環(huán)

callforce(f,en)計(jì)算力

callintegrate(f,en)積分運(yùn)動(dòng)方程

t=t+deltcallsample抽樣平均enddostopend8/16/20245初始化3.1.1初始化要啟動(dòng)模擬，應(yīng)對(duì)體系中所有的粒子賦予初始位置和速度。所選粒子位置應(yīng)與要模擬的結(jié)構(gòu)相容。通常，將粒子初始放在立方晶格上來完成。然后對(duì)每一個(gè)粒子的分速度賦予一個(gè)區(qū)間(-0.5,0.5)上隨機(jī)分布的值。接著，改變所有粒子的速度以確保總動(dòng)量為零，并標(biāo)定速度以使平均動(dòng)能為指定值。熱平衡時(shí)有式中v

是給定粒子速度的

分量。瞬時(shí)溫度T(t):顯然,通過因子fs=[T/T(t)]1/2來調(diào)整所有粒子的速度以使得瞬時(shí)溫度與指定溫度相一致。8/16/20246分子動(dòng)力學(xué)的初始化subroutineinitMD程序的初始化sumv=0sumv2=0doi=1,npart

x(i)=lattice_pos(i)將粒子放在晶格上

v(i)=(ranf()-0.5)給定隨機(jī)速度

sumv=sumv+v(i)質(zhì)心速度

sumv2=sumv2+v(i)*v(i)動(dòng)能enddosumv=sumv/npart

質(zhì)心速度sumv2=sumv2/npart均方速度fs=sqrt(3*temp/sumv2)速度的標(biāo)定因子doi=1,npart

v(i)=(v(i)-sumv)*fs

設(shè)定所需的動(dòng)能和質(zhì)心速度零

xm(i)=x(i)–v(i)*dt

先前時(shí)間步的位置enddoreturnend8/16/20247力的計(jì)算3.1.2力的計(jì)算我們以Lennard-Jones體系為例，來計(jì)算粒子間的互作用力。第i個(gè)原子受的力為：如果粒子對(duì)i，j足夠近以至可發(fā)生作用的話，必須計(jì)算它們之間的力及其對(duì)勢(shì)能的貢獻(xiàn)。假設(shè)希望計(jì)算x方向上的力8/16/20248力的計(jì)算subroutineforce(f,en)確定力和能量en=0doi=1,npart

f(i)=0設(shè)置力為零enddodoi=1,npartdoj=i+1,npart

xr=x(i)-x(j)

xr=xr–box*nint(xr/box)周期性邊界條件

r2=xr*xr

if(r2.lt.rc2)then檢驗(yàn)截?cái)嗑嚯x

r2i=1/r2r6i=r2i**3ff=48*r2i*r6i(r6i-0.5)Lennard-Jones勢(shì)

8/16/20249力的計(jì)算（續(xù)）

f(i)=f(i)+ff*xr

更新力

f(j)=f(j)–ff*xren=en+4*r6i*(r6i-1)-ecut更新能量,ecut=4(rc-12-rc-6)

endif

enddoenddoreturnend8/16/202410運(yùn)動(dòng)方程積分3.1.3運(yùn)動(dòng)方程積分采用Verlet算法，這一算法不僅是最簡(jiǎn)單的，而且通常是最好的。

其推導(dǎo)可通過對(duì)在t時(shí)刻某粒子的坐標(biāo)泰勒展開獲得類似地，將這兩個(gè)方程相加，可得或

(3.1)新位置的估計(jì)值有一t4級(jí)別的誤差，式中t是分子動(dòng)力學(xué)模擬中的時(shí)間步長(zhǎng)。注意Verlet算法并不使用速度來計(jì)算新的位置。8/16/202411運(yùn)動(dòng)方程積分(續(xù))將上面兩個(gè)泰勒展開式相減,可以得到速度(可從軌跡知識(shí)得到)。

或

(3.2)速度的這一表達(dá)式的精度只在t2級(jí)別左右。由上可以看出，如果原子的前兩個(gè)位置是已知的，則根據(jù)式(3.1)和式(3.2)就可進(jìn)行相關(guān)的計(jì)算。然而，在典型的分子動(dòng)力學(xué)模擬中，只有初始位置和初始速度是給定的。所以，在利用遞歸的Verlet算法開始計(jì)算之前，首先要找到原子第二個(gè)位置的計(jì)算方法。一般地，應(yīng)用普通運(yùn)動(dòng)學(xué)方程可得8/16/202412積分動(dòng)動(dòng)方程subroutineintegrate(f,en)積分動(dòng)動(dòng)方程sumv=0sumv2=0doi=1,npartMD循環(huán)

xx=2*x(i)-xm(i)+dt*dt*f(i)/2Verlet算法(3.1)vi=(xx-xm(i))/(2*dt)速度(3.2)

sumv=sumv+vi質(zhì)心速度

sumv2=sumv2+vi*vi

總動(dòng)量

xm(i)=x(i)更新先前步的位置

x(i)=xx更新當(dāng)前步的位置enddotemp=sumv2/(3*npart)瞬時(shí)溫度etot=en+sumv2/2總能量returnend8/16/202413§3.2等壓方法（Isobaricmethod）上節(jié)我們討論的分子動(dòng)力學(xué)是一種研究在體積V中N個(gè)粒子的經(jīng)典體系的自然的時(shí)間演繹的方法。在這些模擬中，總能量E是運(yùn)動(dòng)的常數(shù)。若我們假設(shè)時(shí)間平均等效于系綜平均，那么在一般的MD模擬中得到的時(shí)間平均等價(jià)于在微正則系綜（恒NVE）的系綜。在其他系綜（如NVT或NpT）中進(jìn)行模擬更方便，更接近實(shí)際情況。乍一看，在不是微正則系綜的其他系綜中進(jìn)行MD模擬是不可能。但已提出了兩種相當(dāng)不同的解決這一問題的方法。一種是基于這樣一種思想，即其他系綜的動(dòng)力學(xué)模擬可以通過混合牛頓力學(xué)的MD和一定的MonteCarlo移動(dòng)來實(shí)現(xiàn)。另一種方法完全源自動(dòng)力學(xué)，它基于對(duì)體系拉格朗日運(yùn)動(dòng)方程的重新推導(dǎo)。8/16/202414經(jīng)典粒子服從牛頓力學(xué)、分析力學(xué)？*牛頓力學(xué)可以解決經(jīng)典粒子的運(yùn)動(dòng)問題：牛頓力學(xué)用矢量形式建立運(yùn)動(dòng)基本方程（矢量力學(xué)）。在處理約束系統(tǒng)、變形體的動(dòng)力學(xué)等問題時(shí)有些困難。分析力學(xué)（標(biāo)量力學(xué)）屬于牛頓力學(xué)的范疇，但是又高于牛頓力學(xué)。*分析力學(xué)也可以解決經(jīng)典粒子的運(yùn)動(dòng)問題（哈密頓方程）：qi

是粒子i的廣義坐標(biāo)，pi是與之共軛的廣義動(dòng)量。8/16/202415牛頓力學(xué)用力、加速度等矢量表示質(zhì)點(diǎn)運(yùn)動(dòng)規(guī)律，著眼于體系內(nèi)每個(gè)質(zhì)點(diǎn)受到的力及其產(chǎn)生的效果。分析力學(xué)通過動(dòng)能、勢(shì)能等標(biāo)量表示運(yùn)動(dòng)規(guī)律，并能將著眼點(diǎn)由質(zhì)點(diǎn)轉(zhuǎn)為體系整體（能量是標(biāo)量，是