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

1、引言分子動(dòng)力學(xué)方法(Molecular計(jì)算機(jī)模擬中的另一類確定性模擬方法，即統(tǒng)計(jì)物理中的所謂合于動(dòng)力學(xué)方法Dynamics Method ）。這種方法是按該體系內(nèi)部的內(nèi)稟動(dòng)力學(xué)規(guī)律來計(jì)算并確定位形的轉(zhuǎn)變。它首先需要得到每個(gè)時(shí)刻各建立一組分子的運(yùn)動(dòng)方程，并通過直接對(duì)系統(tǒng)中的一個(gè)個(gè)分子運(yùn)動(dòng)方程進(jìn)行數(shù)值求解, 個(gè)分子的坐標(biāo)與動(dòng)量，即在相空間的運(yùn)動(dòng)軌跡，再利用統(tǒng)計(jì)計(jì)算方法得到多體系統(tǒng)的靜態(tài)和動(dòng)態(tài)特性,MD方法中不存在任何隨機(jī)因素。在從而得到系統(tǒng)的宏觀性質(zhì)。在這樣的處理過程中我們可以看出:MD方法處理過程中方程組的建立是通過對(duì)物理體系的微觀數(shù)學(xué)描述給出的。在這個(gè)微觀的物理體系 中，每個(gè)分子都各自服從經(jīng)典

2、的牛頓力學(xué)。每個(gè)分子運(yùn)動(dòng)的內(nèi)稟動(dòng)力學(xué)是用理論力學(xué)上的哈密頓量或者拉格朗日量來描述，也可以直接用牛頓運(yùn)動(dòng)方程來描述。確定性方法是實(shí)現(xiàn) Boltzman的統(tǒng)計(jì)力學(xué)途徑。這種方法可以處理與時(shí)間有關(guān)的過程，因而可以處理非平衡態(tài)問題。但是使用該方法的程序較復(fù)雜，討算量大，占內(nèi)存也多、本節(jié)將介紹分子動(dòng)力學(xué)方法及其應(yīng)用。原則上，MD方法所適用的微觀物理體系并無什么限制。這個(gè)方法適用的體系既可以是少體系統(tǒng)，也可以是多體系統(tǒng)；既可以是點(diǎn)粒子體系，也可以是具有內(nèi)部結(jié)構(gòu)的體系；處理的微觀客體既可以是分子，也可以是其他的微觀粒子。實(shí)際上，MD模擬方法和隨機(jī)模擬方法一樣都面臨著兩個(gè)基本限制：一個(gè)是有限觀測(cè)時(shí)間的限制；另

3、一個(gè)是有限系統(tǒng)大小的限制。通常人們感興趣的是體系在熱力學(xué)極限下（即粒子數(shù)日趨于無窮時(shí)）的性質(zhì)。但是計(jì)算機(jī)模擬允許的體系大小要比熱力學(xué)極限小得多，因此可能會(huì)出現(xiàn)有限尺寸效應(yīng)。為了減小有限尺寸效應(yīng)，人們往往引入周期性、全反射、漫反射等邊界條件。 當(dāng)然邊界條件的引入顯然會(huì)影響體系的某些性質(zhì)。對(duì)于MD方法，向然的系綜是微正則系綜，這時(shí)能量是運(yùn)動(dòng)常量。然而，當(dāng)我們想要研究溫度和（或）壓力是運(yùn)動(dòng)常量的系統(tǒng)時(shí)，系統(tǒng)不再是封閉的。例如當(dāng)溫度為常量的系統(tǒng)可以認(rèn)為系統(tǒng)是放置在 一個(gè)熱俗中。當(dāng)然，在 MD方法中我們只是在想像中將系統(tǒng)放入熱浴中。實(shí)際上，在模擬計(jì)算中具體 所采取的做法是對(duì)一些自由度加以約束。例如在恒溫

4、體系的情況下， 體系的平均動(dòng)能是一個(gè)不變量。這時(shí)我們可以設(shè)計(jì)一個(gè)算法， 使平均動(dòng)能被約束在一個(gè)給定值上。由于這個(gè)約束，我們并不是在真正處理一個(gè)正則系綜，而實(shí)際上僅僅是復(fù)制了這個(gè)系綜的位形部分。只要這一約束不破壞從一個(gè)狀態(tài)到另一個(gè)狀態(tài)的馬爾科夫特性，這種做法就是正確的。不過其動(dòng)力學(xué)性質(zhì)可能會(huì)受到這一約束的影響。自五十年代中期開始，MD方法得到了廣泛的應(yīng)用。 它與蒙特卡洛方法一起已經(jīng)成為計(jì)算機(jī)模擬的重要方法。應(yīng)用 MD方法取得了許多重要成果，例如氣體或液體的狀態(tài)方程、相變問題、吸附問題等， 以及非平衡過程的研究。其應(yīng)用已從化學(xué)反應(yīng)、生物學(xué)的蛋白質(zhì)到重離子碰撞等廣泛的學(xué)科研究領(lǐng)域。、分子運(yùn)動(dòng)方程的數(shù)

5、值求解采用MD方法時(shí)，必須對(duì)一組分于運(yùn)動(dòng)微分方程做數(shù)值求解。從計(jì)算數(shù)學(xué)的角度來看，這個(gè)求解是一個(gè)初值問題。實(shí)際上計(jì)算數(shù)學(xué)為了求解這種問題己經(jīng)發(fā)展了許多的算法，但并不是所有的這些算法都可以用來解決物理問題。下面我們先以一個(gè)一維諧振子為例，來看一下如何用計(jì)算機(jī)數(shù)值計(jì)算方法求 解初值問題。一維諧振子的經(jīng)典哈密頓量為(2.1)這里的哈密頓量（即能量）為守恒量。假定初始條件為x（p）、p（0）,則它的哈密頓方程是對(duì)時(shí)間的一階微分方程<Lr_ pdim凸旅（2.2）現(xiàn)在我們要用數(shù)值積分方法計(jì)算在相空間中的運(yùn)動(dòng)軌跡（X（t）、p（t）。我們采用有限差分法，將微分方程變?yōu)橛邢薏罘址匠?，以便在?jì)算機(jī)上做數(shù)

6、值求解，并得到空間坐標(biāo)和動(dòng)量隨時(shí)間的演化關(guān)系。首先，我們?nèi)〔罘钟?jì)算的時(shí)間步長為h，采用我們有限差分法中的一階微分形式的向前差商表示，即直接運(yùn)用展開到h的一階泰勒展開公式即(2.3)則微分方程(2.2)可以被改寫為差分形(2.4)(2.5)將上面兩個(gè)公式整理后，我們得到解微分方程(2.2)的歐拉(Euler)算法如砂5+磅錯(cuò)存咖)(2.9)帝)(2.10)(2.6)加 + 防二 p(t)-hkx(t)( 2.7)這是x(t)、p(t)的一組遞推公式、有了初始條件 x(0)、p(0)，就可以一步一步地使用前一時(shí)刻的坐標(biāo)、動(dòng) 量值確定下一時(shí)刻的坐標(biāo)、動(dòng)量值。這個(gè)方法是一步法的典型例子。由于在實(shí)際數(shù)值

7、計(jì)算時(shí) h的大小是有限的，因而在上述算法中微分被離散化為差分形式來計(jì)算時(shí)總 是有誤差的?？梢宰C明一步法的局部離散化誤差與總體誤差是相等的，都為0(h2)的量級(jí)。在實(shí)際應(yīng)用中，適當(dāng)?shù)剡x擇h的大小是十分重要的。h取得太大，得到的結(jié)果偏離也大，甚至于連能量都不守恒：h取得太小，有可能結(jié)果仍然不夠好。這就要求我們改進(jìn)計(jì)算方法，進(jìn)一步考慮二步法。實(shí)際上泰勒展開式的一般形式代+山/(0+ f4- 0(沖)訊 1-(2.8)其中0(hn+1)表示誤差的數(shù)量級(jí)。前面敘述的歐拉算法就是取n=1?，F(xiàn)在考慮公式(2.8)中直到含h的二次項(xiàng)的展開(即取 n= 2),則得到將上面兩式相加、減得到含二階和一階導(dǎo)數(shù)的公式乎

8、=+廿 + 4 2/(o+f(f -曲(2.11)(2.12)令f(t) = x(t)，利用牛頓第二定律公式F(t)吩，公式(2.11)寫為坐標(biāo)的遞推公式W十利二-X(t -h) + 2式)十護(hù)型"(2.13)公式(2.12)寫為計(jì)算動(dòng)量的公式得到P(r)=就 G) = mv(i) = r(r + A)-升)(2.14)這樣我們就推導(dǎo)出了一個(gè)比(2.6)和(2.7)更精確的遞推公式。這是二步法的一種，稱為Verle方法。還有其他一些二步法，如龍格一庫塔(Ru nge Kutta )方法等。當(dāng)然我們還可以建立更高階的多步算法，然而大部分更高階的方法所需要的內(nèi)存比一步法和二步法所需要的大

9、得多，并且有些更高階的方法還需要用迭代來解出隱式給定的變量，內(nèi)存的需求量就更大。 并且當(dāng)今的計(jì)算機(jī)都僅僅只有有限的內(nèi)存，因而并不是所有的高階算法都適用于物理系統(tǒng)的計(jì)算機(jī)計(jì) 算。在實(shí)際數(shù)值計(jì)算中，我們必須特別注意舍入誤差和穩(wěn)定性問題。為了減少舍入誤差，我們可以采用高精度計(jì)算，并且要避免相近大小的數(shù)相消，以及數(shù)量級(jí)相差很大的兩個(gè)數(shù)相加和注意運(yùn)算順序。三、分子動(dòng)力學(xué)模擬的基本步驟在計(jì)算機(jī)上對(duì)分子系統(tǒng)的 MD模擬的實(shí)際步驟可以劃分為四步；首先是設(shè)定模擬所采用的模型: 第二，給定初始條件；第三，趨于平衡的計(jì)算過程；最后是宏觀物理量的計(jì)算。下面就這四個(gè)步驟分別 做簡(jiǎn)單介紹。1、模擬模型的設(shè)定設(shè)定模型是分子

10、動(dòng)力學(xué)模擬的第一步工作。例如在一個(gè)分子系統(tǒng)中， 假定兩個(gè)分子間的相互作用勢(shì)為硬球勢(shì)，其勢(shì)函數(shù)表示為如果如果實(shí)際上，更常用的是圖(3.1) Lennard-Jones型勢(shì)。它的勢(shì)函數(shù)表示為V(r)=4e(3.1)其中-£是位勢(shì)的最小值 確定為長度的單位)(&可以確定能量的單位),這個(gè)最小值出現(xiàn)在距離r等于21/6c的地方(c可以0.6 0.8 1. 0 1.2 t.4 M E 2.0 2.2 2r4 2.6圖 3.1 Lennard-Jones 勢(shì)例如對(duì)在微正則系綜的模模型確定后，根據(jù)經(jīng)典物理學(xué)的規(guī)律我們就可以知道在系綜模擬中的守恒量。 擬中能量、動(dòng)量和角動(dòng)量均為守恒量。在此系

11、綜中它們分別表示為(3.2)i（3.3）i（3.4）其中Pi=m門。由于我們只限于研究大塊物質(zhì)在給定密度卜的性質(zhì)，所以必須引進(jìn)一個(gè)叫做分子動(dòng)力學(xué)元胞的體積元，以維持一個(gè)恒定的密度。對(duì)氣體和液體，如果所占體積足夠大，并且系統(tǒng)處于熱平衡狀態(tài) 的情況下，那么這個(gè)體積的形狀是無關(guān)緊要的。對(duì)于晶態(tài)的系統(tǒng)，元胞的形狀是有影響的。 為了計(jì)算簡(jiǎn)便，對(duì)于氣體和液體，我們?nèi)∫粋€(gè)立方形的體積為MD元胞。設(shè)MD元胞的線度大小為 L,則其體積為L3。由于引進(jìn)這樣的立方體箱子，將產(chǎn)生六個(gè)我們不希望出現(xiàn)的表面。模擬中碰撞這些箱的表面的 粒子應(yīng)當(dāng)被反射回到元胞內(nèi)部，特別是對(duì)粒子數(shù)目很少的系統(tǒng)。然而這些表面的存在對(duì)系統(tǒng)的任何一

12、種性質(zhì)都會(huì)有重大的影響。 為了減小引入的表面效應(yīng)，我們采用周期性邊界條件。采用這種邊界條件，我們就可以消除引入的表面效應(yīng)，構(gòu)造出一個(gè)準(zhǔn)無窮大的體積來更精確地代表宏觀系統(tǒng)。實(shí)際上，這里我們做了一個(gè)假定，即讓這個(gè)小體積元胞鑲嵌在一個(gè)無窮大的大塊物質(zhì)之中。周期性邊界條件的數(shù)學(xué)表示形式為A（x） - A（xnL）t 2=（%代厲）（3.5）其中A為任意的可觀測(cè)量，nn2、n3。為任意整數(shù)。這個(gè)邊界條件就是命令基本 MD元胞完全等同地 重復(fù)無窮多次。具體在實(shí)現(xiàn)該邊界條件時(shí)是這樣操作的：當(dāng)有一個(gè)粒子穿過基本MD元胞的六方體表面時(shí)，就讓這個(gè)粒子以相同的速度穿過此表面對(duì)面的表面重新進(jìn)入該MD元胞內(nèi)。在分子動(dòng)力

13、學(xué)模擬中考慮粒子間的相互作用時(shí)，通常采用最小像力約定。這個(gè)約定是在由無窮重復(fù)的MD基本元胞中，一個(gè)粒子只同它所在的基本元胞內(nèi)的另外N-1個(gè)（設(shè)在此元胞內(nèi)有 N個(gè)粒子）中的每個(gè)粒子或其最鄰近影像粒子發(fā)生相互作用。如果ri處的粒子i同rj處的粒子j之間的距離為（對(duì)一切的 n）（3.6）實(shí)際上這個(gè)約定就是通過滿足不等式條件rc<L/2來截?cái)辔粍?shì)（L為截止距離）。通常L的數(shù)值應(yīng)當(dāng)選得很大，使得距離大于 L/2的粒子的相互作用可以忽略，以避免有限尺寸效應(yīng)。采用最小像力約定使得在截?cái)嗵幜Ｗ拥氖芰τ幸粋€(gè)函數(shù)的奇異性，這會(huì)給模擬計(jì)算帶來誤差。2、給定初始條件MD模擬進(jìn)人對(duì)系統(tǒng)微分方程組做數(shù)值求解的過程時(shí)

14、，需要知道粒子的初始位置和速度的數(shù)值。不同的算法要求不同的初始條件。例如，Verlet方法需要兩組坐標(biāo)來啟動(dòng)計(jì)算；一組是零時(shí)刻的坐標(biāo)，另一組是前進(jìn)一個(gè)時(shí)間步長時(shí)的坐標(biāo)，或者是一組零時(shí)刻的速度值。但是，一般來說系統(tǒng)的初始條件都是不可能知道的。表面上看這是一個(gè)難題。實(shí)際上，精確選擇待求系統(tǒng)的初始條件是沒有什么意義的，因?yàn)槟M時(shí)間足夠長時(shí)，系統(tǒng)就會(huì)忘掉初始條件。但是初始條件的合理選擇將可以加快系統(tǒng)趨于平衡。常用的初始條件可以選擇為：（1）令初始位置在差分劃分網(wǎng)格的格子上，初始速度則從玻爾茲曼分布隨機(jī)抽樣得到。（2 ）令初始位置隨機(jī)地偏離差分劃分網(wǎng)格的格子，初始速度為零。（3）令初始位置隨機(jī)地偏離差分

15、劃分網(wǎng)格的格子，初始速度從玻爾茲曼分布隨機(jī)抽樣得到。3趨于平衡按照上面給出的運(yùn)動(dòng)方程、邊界條件和初始條件，就可以進(jìn)行分子動(dòng)力學(xué)模擬計(jì)算。但是，這樣計(jì)算出的系統(tǒng)不會(huì)具有所要求的系統(tǒng)能量，并且這個(gè)狀態(tài)本身也還不是一個(gè)平衡態(tài)。為了使系統(tǒng)達(dá)到平衡，模擬中需要一個(gè)趨衡過程。在這個(gè)過程中，增加或從系統(tǒng)中移出能量，直到系統(tǒng)具有所要求的能量。然后，再對(duì)運(yùn)動(dòng)方程中的時(shí)間向前積分若于步，使系統(tǒng)持續(xù)給出確定能量值。我們稱；這時(shí)系統(tǒng)已經(jīng)達(dá)到平衡態(tài)。這段達(dá)到平衡所需的時(shí)間稱為弛豫時(shí)間。在MD模擬中，時(shí)間步長 h的大小選擇是十分重要的。它決定了模擬所需要的時(shí)間。為了減小誤差，步長h必須取得小一些；但是取得太小，系統(tǒng)模擬的

16、弛豫時(shí)間就很長。這里需要積累一定的模擬經(jīng)驗(yàn)，選擇適當(dāng)?shù)臅r(shí)間步長h。例如，對(duì)一個(gè)具有幾百個(gè)氬（Ar）分于的體系，如果采用Lennard-Jones位勢(shì)，我們發(fā)現(xiàn)取h為10-2量級(jí)，就可以得到好的相圖這里 選擇的h是沒有量綱的，實(shí)際上這樣選擇的h對(duì)應(yīng)的時(shí)間在10-14秒的量級(jí)。如果模擬1000步，系統(tǒng)達(dá)到平衡態(tài)，馳豫時(shí)間只有 10-11秒。4宏觀物理量的計(jì)算實(shí)際計(jì)算宏觀物理量往往是在MD模擬的最后階段進(jìn)行的。它是沿著相空間軌跡求平均來計(jì)算得到的。例如對(duì)于一個(gè)宏觀物理量A，它的測(cè)量值應(yīng)當(dāng)為平均值A(chǔ)。如果已知初始位置和動(dòng)量為 r(N)(o)和P(N)(0)(上標(biāo)N表示系綜N個(gè)粒子的對(duì)應(yīng)坐標(biāo)和動(dòng)量參數(shù))

17、，選擇某種MD算法求解具有初值問題 的運(yùn)動(dòng)方程，便得到相空間軌跡 r(N)(t)， p(N)(t)對(duì)軌跡平均的宏觀物理量A的表示為兀阿f切(" 心林約)t卄社一 4八(3.7)如果宏觀物理量為動(dòng)能，它的平均為取二您(妣匕(#“&)3.8)可以表示為在時(shí)間的各(3.8)由于在模擬過程中計(jì)算出的動(dòng)能值是在不連續(xù)的路徑上的值，因此公式(個(gè)間斷點(diǎn)卩上計(jì)算動(dòng)能的平均值(3.9)在MD模擬過程中，溫度是需要加以監(jiān)測(cè)的物理量, 我們可以從平均動(dòng)能值計(jì)算得到溫度值：特別是在模擬的起始階段。根據(jù)能量均分定理,(3.12)(3.10)其中d為每個(gè)粒子的自由度，如果不考慮系統(tǒng)所受的約束，則d= 3

18、。系統(tǒng)內(nèi)部的位形能量的軌道平均值為(3.11)nnQ尸心假定位勢(shì)在rc處被截?cái)?，那么上式?jì)算出的勢(shì)能以及由此得到的總能量就包含有誤差。為了對(duì)此偏差 作出修正，我們采用對(duì)關(guān)聯(lián)函數(shù)來表示位能式中的g(r)就是對(duì)關(guān)聯(lián)函數(shù)，它是描述與時(shí)間無關(guān)的粒子間關(guān)聯(lián)性的量度。g(r)的物理意義是當(dāng)原點(diǎn)r=0處有一個(gè)粒子時(shí)，在空間位置r的點(diǎn)周圍的體積元中單位體積內(nèi)發(fā)現(xiàn)另一個(gè)粒子的幾率。若n(r)為距離原點(diǎn)r到r+ Ar之間的平均粒子數(shù)，則V n(r)(3.13)N在MD模擬過程中，所有的距離已經(jīng)在力的計(jì)算中得到，因而很容易計(jì)算對(duì)關(guān)聯(lián)函數(shù)的值。圖3.2為由計(jì)算機(jī)模擬得到的兩組不同參數(shù)下的對(duì)關(guān)聯(lián)函數(shù)的例子。由于位勢(shì)的截

19、斷，對(duì)關(guān)聯(lián)函數(shù)僅對(duì)rc< L/2以下的距離有意義。在公式(3.11 )中，所有的位能都加到截?cái)嗑嚯x為止，尾部修正可以取為Uc - 2np j «(r)g(r)r2dr(3.14)壓強(qiáng)可以通過計(jì)算在面積元 dA的法線方向上凈動(dòng)量轉(zhuǎn)移的時(shí)間平均值來得到，也可以利用含對(duì)關(guān)聯(lián)函 數(shù)的維里狀態(tài)方程計(jì)算。該維里狀態(tài)方程可以寫為p =旳TJ7 S（r）字4”擊（3.15）一項(xiàng)是對(duì)位至于勢(shì)能的計(jì)算，我們可以把積分劃分為兩項(xiàng)，一項(xiàng)是由相互作用力程之內(nèi)的貢獻(xiàn)引起的,勢(shì)截?cái)嗟母恼?xiàng)：(3.16)其中長程改正項(xiàng)為:二牛 J7 g（廣）尋 4r3dr(3.17)F面我們將討論具體如何進(jìn)行 MD模擬:r/

20、a圖6.2由計(jì)算機(jī)模擬得到的兩組不同參數(shù)下的對(duì)關(guān)聯(lián)函數(shù)4、平衡態(tài)分子動(dòng)力學(xué)模擬在經(jīng)典MD模擬方法的應(yīng)用當(dāng)中，存在著對(duì)兩種系統(tǒng)狀態(tài)的MD模擬。一種是對(duì)平衡態(tài)的 MD模擬，另一類是對(duì)非平衡態(tài)的 MD模擬。對(duì)平衡態(tài)系綜 MD模擬又可以分為如下類型： 微正則系綜的 MD （NVE）模擬，正則系綜的MD （NVT ）模擬，等溫等壓系綜 MD（ NPT）模擬和等焓等壓系綜 MD （NPH） 模擬等。下面我們僅對(duì)平衡態(tài)的MD方法中前兩類模擬做簡(jiǎn)單的介紹。1）、微正則系綜的 MD模擬在進(jìn)行對(duì)微正則系綜的MD模擬時(shí)，首先我們要確定所采用的相互作用模型。我們假定一個(gè)孤立的多粒子體系，其粒子間的相互作用位勢(shì)是球?qū)ΨQ

21、的，則其哈密頓量可以寫為(4.1其中rij是第i個(gè)粒子與第j個(gè)粒子之間的距離。在這個(gè)微正則系綜中，由于這個(gè)系統(tǒng)的哈密頓量中不顯式地出現(xiàn)時(shí)間關(guān)聯(lián)，因而系統(tǒng)的能量是個(gè)守恒量。系統(tǒng)的體積和粒子數(shù)也是不變的。此外；由于整個(gè)系統(tǒng)并未運(yùn)動(dòng)，所以整個(gè)系統(tǒng)的總動(dòng)量P恒等于零。這就是系統(tǒng)受到的四個(gè)約束。由該系統(tǒng)的哈密頓量可以推導(dǎo)出牛頓方程形式的運(yùn)動(dòng)方程組dtm疋,N)(4.2)要用數(shù)值求解的方法解出（ 求解變成求解方程組：4.2）微分方程組，類似于本章第二節(jié)中介紹的Veriet方法方程組（4.2）的片(f + A) - 2rr(z) - r/r - A) +/刑，(d =丄2卄,N)該方程組反映出：從前面t和t

22、-h時(shí)刻這兩步的空間坐標(biāo)位置及t時(shí)刻的作用力，就可以算出下一步 t+h時(shí)刻的坐標(biāo)位置。下面為了將(4.3)式寫成更簡(jiǎn)潔的形式，令(4.4)則從(4.3)式可以得到如下差分方程組的形式嚴(yán) 1)="嚴(yán)一嚴(yán) + f %/mt“ 中(4.5)如果已知一組初始空間位置ri(0),ri(1)，則通過求解方程組(4.5) 步步得到r,2) ,ri(3),。由空間坐標(biāo)又可以算出粒子的運(yùn)動(dòng)速度為(4.6)這里在第n+ I步算出的速度是前一時(shí)刻； 即第n步的速度、因而動(dòng)能的計(jì)算比勢(shì)能的計(jì)算落后一步。根據(jù)上述原理我們可以將粒子數(shù)恒定、體積恒定、能量恒定的微正則系綜(NVE )的MD模擬步驟設(shè)計(jì)如下：(1)

23、 給定初始空間位置 $0), ) , (i=1,2,3,N)(2) 計(jì)算在第n步時(shí)粒子所受的力Fi(n) : Fi嘰 Fg)。(3) 利用公式ri(n12ri(nri(nJ)Fi(n)h2 /m計(jì)算在第n+l步時(shí)所有粒子所處的空間位置ri(n1)(4) 計(jì)算第 n 步的速度：Vi(n) =(ri(n1)-ri(nJ)/2h(5) 返回到步驟(2),開始下一步的模擬計(jì)算。如前所述，用上述形式的Verlet算法，動(dòng)能的計(jì)算比勢(shì)能的計(jì)算落后一步。此外，這種算法不是自啟動(dòng)的。要真正求出微分方程組(4.2)的解，除了需要給出初始空間位置ri(0)夕卜，還要求另外給出一組空間位置ri(1)。實(shí)際，有時(shí)候采

24、用改進(jìn)后的計(jì)算方法可能更方便；即把N個(gè)粒于的初始位置放在網(wǎng)格的格點(diǎn)上，然后加以擾動(dòng)。如果初始條件是空間位置和速度，則采用下面的公式來計(jì)算空間位置n(1)(4.7)然后再按上述模擬步驟進(jìn)行計(jì)算。verlet算法的速度變型形式將會(huì)使其數(shù)值計(jì)算的穩(wěn)定性得到加強(qiáng)。下面我們就此做簡(jiǎn)單介紹。令則公式(4.5)寫為(4.9)Verlet算法的速度形式的模擬步驟可以表上式在數(shù)學(xué)上與(4.5)式是等價(jià)的，并稱為相加形式。由此 述為(I)給定初始空間位置ri(1) , (i=1,2,3,N)(2) 給定初始速度Vi(1)(3) 利用公式：r," 冷嚴(yán)-hv(n) - Fi(n)h2/2m，計(jì)算在第n +

25、 1步時(shí)所有粒子所處的空間位置(4)計(jì)算在第 n+1 步時(shí)所有粒子的速度v(n1) :v(n 1) =v(n) - h(Fi(n1) - Fi(n)/2m(5)返回到步驟(3)，開始第n+2步的模擬計(jì)算。Verlet速度形式的算法比前一種算法好些。它不僅可以在計(jì)算中得到同一時(shí)間步上的空間位置和速 度，井且數(shù)值計(jì)算的穩(wěn)定性也提高了。一般情況下，對(duì)于給定能量的系統(tǒng)不可能給出精確的初始條件。這時(shí)需要先給出一個(gè)合理的初始條件，然后在模擬過程中逐漸調(diào)節(jié)系統(tǒng)能量達(dá)到給定值。其步驟為：首先將運(yùn)動(dòng)方程組解出若干步的結(jié)果；然后計(jì)算出動(dòng)能和位能；假如總能量不等于給定恒定值，則通過對(duì)速度的調(diào)整來實(shí)現(xiàn)能量守恒。也就是

26、然后再回到第一步，對(duì)下一時(shí)刻的運(yùn)動(dòng)方程求解。反復(fù)進(jìn)行上面的過程，直到系統(tǒng)達(dá)到平衡。這樣的模擬過程也稱為平衡化階段。采用對(duì)速度標(biāo)度的辦法， 可以使速度發(fā)生很大變化。 為了消除可能帶來的效應(yīng)， 必須要有足夠的時(shí) 間讓系統(tǒng)再次建立平衡。在到達(dá)趨衡階段以后，必須檢驗(yàn)粒子的速度分布是否符合麥克斯韋-波爾茲曼分布。2)、正則系綜的MD模擬在統(tǒng)計(jì)物理中的正則系綜模擬是針對(duì)一個(gè)粒子數(shù)N、體積V、溫度T和總動(dòng)量(Pi =0 )i為守恒量的系綜(NVT )。這種情況就如同一個(gè)系統(tǒng)置于熱浴之中，此時(shí)系統(tǒng)的能量可能有漲落，但系 統(tǒng)溫度則已經(jīng)保持恒定。 在正則系綜的 MD模擬中施加的約束與微正則系綜中的不一樣。正則系綜

27、MD方法是在運(yùn)動(dòng)方程組上加上動(dòng)能恒定(即溫度恒定)的約束，而不是像做正則系綜的MD模擬中對(duì)運(yùn)動(dòng)方程加上能量恒定的約束。在正則系綜MD的平衡化過程中，速度標(biāo)度因子一般選下面的形式較為我們可將正則系綜 MD的Verlet算法的速度形式的模擬具體步驟列在下面:(I)給定初始空間位置t(1)：(i =1,2,N)(2) 給定初始速度Vi(3) 利用公式r/n*-hv(n) Fi(n)h2/2m計(jì)算在第n+1步時(shí)所有粒子所處的空間位置ri(n1)(4)計(jì)算在第n + 1步時(shí)所有粒于的速度v(n 1) =v(n) - h(Fi(n 1) - Fi(n)/2m，動(dòng)能和速度標(biāo)度因子| (3N - +)iT.i

28、n+1步粒子的速度(5)計(jì)算將速度 Vin 1標(biāo)度因于 3值，并讓該值作為下一次計(jì)算時(shí)，第廠 >Wn1:(6)返回到步驟(3),開始第n + 2步的模擬計(jì)算。按照上面的步驟，對(duì)時(shí)間進(jìn)行一步步的循環(huán)。待系統(tǒng)達(dá)到平衡后， 則退出循環(huán)。這就是正則系綜的MD模擬過程。下面我們舉一個(gè)微正則系綜的 MD模擬的應(yīng)用示例來看看模擬的結(jié)果。 例： 對(duì)一個(gè)總能量確定的單原于(氬)粒子系統(tǒng)的MD模擬計(jì)算。我們具體選取256個(gè)原子的模擬。粒于間的相互作用位勢(shì)為Lennard-Jones勢(shì)(4.12)其中-£為位勢(shì)的極小值(取&為能量單位)，其位置在r=21/6 b處、該體系的粒于限制在一個(gè)立方體的箱子，邊界上采用最小象力約定。我們采用自然單位制。長度和時(shí)間的標(biāo)度單位分別為b和(m b/48 &1/2(對(duì)氬原于該時(shí)間單位為3 X10-12秒),這樣就使得運(yùn)動(dòng)方程為無量綱形式。模擬時(shí)我們考慮兩個(gè)相圖上的點(diǎn)(T* , p*)=(2.53 , 0.636), ( 0.722 , 0.83134),它們分別具有兩種立方體的尺寸，即L = 7.83和L = 6.75。初始條件假定為：各個(gè)原子處于一個(gè)面心立方格子的格點(diǎn)上，而速度按相應(yīng)溫度下的波爾茲曼分布抽樣取值。位勢(shì)的截?cái)嗳?/p>