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

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

2、系中， 每個(gè)分子都各自服從經(jīng)典的牛頓力學(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ī)模擬方

3、法一樣都面臨著兩個(gè)基本限制： 一個(gè)是有限觀測(cè)時(shí)間的 限制； 另一個(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 方法中我們只是在想像中將

4、系統(tǒng)放入熱浴中。 實(shí)際上， 在模擬計(jì)算中具體所采取的做法是對(duì)一些自由度加以約束。例如在恒溫體系的情況下， 體系的平均動(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)方程、 相變問題

5、、 吸附問題等， 以及非平衡過程的研究。 其應(yīng)用已從化學(xué)反應(yīng)、 生物學(xué)的蛋白質(zhì)到重離子碰撞等廣 泛的學(xué)科研究領(lǐng)域。二、分子運(yùn)動(dòng)方程的數(shù)值求解采用 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(2.1)x(p)、p(0),則它的哈密頓方程是對(duì)時(shí)間的這里的哈密頓量(即能量)為守恒量。假定初始條件為 一階微分方程dx BH _ p

6、d/(2.2)我們采用有限差分法，將現(xiàn)在我們要用數(shù)值積分方法計(jì)算在相空間中的運(yùn)動(dòng)軌跡(X(t)、p(t)微分方程變?yōu)橛邢薏罘址匠?，以便在?jì)算機(jī)上做數(shù)值求解，并得到空間坐標(biāo)和動(dòng)量隨時(shí)間的演化關(guān)系。首先，1我們?nèi)〔罘钟?jì)算的時(shí)間步長(zhǎng)為h,采用我們有限差分法中的一階微分形式的向前差商表示，即直接運(yùn)用展開到h的一階泰勒展開公式+ /i) = /|t) + h + 0(用)<y二九必卜并) drh(2.3)<ix 過+血)_呦_加)則微分方程(2.2)可以被改寫為差分形(2.4)dr(2.5)將上面兩個(gè)公式整理后，我們得到解微分方程(2.2)的歐拉(Euler)算法W + 由)=x(z)+ 如

7、(x)m(2.6) 衛(wèi)C十人)二(2.7)這是x(t)、p(t)的一組遞推公式、有了初始條件 x(0)、p(0)，就可以一步一步地使用前一時(shí)刻的坐 標(biāo)、動(dòng)量值確定下一時(shí)刻的坐標(biāo)、動(dòng)量值。這個(gè)方法是一步法的典型例子。由于在實(shí)際數(shù)值計(jì)算時(shí) h的大小是有限的，因而在上述算法中微分被離散化為差分形式來計(jì) 算時(shí)總2)的量級(jí)。在實(shí)際應(yīng)用是有誤差的。 可以證明一步法的局部離散化誤差與總體誤差是相等 的，都為0(h中，適當(dāng)?shù)剡x擇 h的大小是十分重要的。h取得太大，得到的結(jié)果偏離也大，甚至 于連能量都不守恒：h取得太小，有可能結(jié)果仍然不夠好。這就要求我們改進(jìn)計(jì)算方法，進(jìn)一步 考慮二步法。實(shí)際上泰勒展開式的一般形

8、式n 石I+ 0的)(2.8)n+1)表示誤差的數(shù)量級(jí)。前面敘述的歐拉算法就是取n=1?，F(xiàn)在考慮公式(2.8)中直到含h其中O(h的二次項(xiàng)的展開(即取 n = 2)，則得到d? 2 dr f(2.9)/(t - a) = /(/)-+ - + o3)d/ 2 dx2 ' '(2.7) 將上面兩式相加、減得到含二階和一階導(dǎo)數(shù)的公式dz為硬球勢(shì)，其勢(shì)函數(shù)表示為L(zhǎng)enn ard-Jo nes型勢(shì)。它的勢(shì)函數(shù)表示為實(shí)際上，更常用的是圖(3.1)2h(2.11)(2.12)2xdF(t) m,公式，利用牛頓第二定律公式 x(t)f(t)令=(2.11)寫為坐標(biāo)的遞推公式 2dt 2十血)

9、=一* 一 ) + 2加)+護(hù)旦。(2.13)(2.12) 寫為計(jì)算動(dòng)量的公式得到公式Zh血)=/nx(r) = mv(f)=曲 + /r) - x(f - A)(2.14)方法。還有更精確的遞推公式。這是二步法的一種，稱為Verle這樣我們就推導(dǎo)岀了一個(gè)比 (2.6)和(2.7) Kutta )方法等。其他一些二步法，如龍格一庫(kù)塔(Runge然而大部分更高階的方法所需要的內(nèi)存比一步法和二步法當(dāng)然我們還可以建立更高階的多步算法，所需要的大得多，并且有些更高階的方法還需要用迭代來解岀隱式給定的變量，內(nèi)存的需求量就更大。并且當(dāng)今的計(jì)算機(jī)都僅僅只有有限的內(nèi)存，因而并不是所有的高階算法都適用于物理系統(tǒng)

10、的計(jì)算機(jī)計(jì)算。我們可以采用為了減少舍入誤差，在實(shí)際數(shù)值計(jì)算中，我們必須特別注意舍入誤差和穩(wěn)定性問題。并且要避免相近大小的數(shù)相消，以及數(shù)量級(jí)相差很大的兩個(gè)數(shù)相加和注意運(yùn)算順序。高精度計(jì)算，三、 分子動(dòng)力學(xué)模擬的基本步驟 模擬的實(shí)際步驟可以劃分為四步；首先是設(shè)定模擬所采用的模型：MD 在計(jì)算機(jī)上對(duì)分子系統(tǒng)的第二，給定初始條件；第三，趨于平衡的計(jì)算過程；最后是宏觀物理量的計(jì)算。下面就這四個(gè)步 驟分別做簡(jiǎn)單介紹。、模擬模型的設(shè)定1假定兩個(gè)分子間的相互作用勢(shì)設(shè)定模型是分子動(dòng)力學(xué)模擬的第一步工作。例如在一個(gè)分子系統(tǒng)中，V(r) = 4f<yr）（3.11/6可以b的地方b（ 2,（- £是

11、位勢(shì)的最小值£可以確定能量的單位）這個(gè)最小值岀現(xiàn)在距離等于其中確定為長(zhǎng)度的單位）3.1 Lennard-Jo nes勢(shì)圖例如對(duì)在微正則系綜的模根據(jù)經(jīng)典物理學(xué)的規(guī)律我們就可以知道在系綜 模擬中的守恒量。模型確定后，擬中能量、動(dòng)量和角動(dòng)量均為守恒量。在此系綜中它們分別表護(hù)+也）示為)(3.2 3(3.33.4）（所以必須引進(jìn)一個(gè)叫做分子動(dòng)力學(xué)元。由于我們只限于研究大塊物質(zhì)在給定密度卜的性質(zhì)，其中p=mrii胞的體積元，以維持一個(gè)恒定的密度。對(duì)氣體和液體，如果所占體積足夠大，并且系統(tǒng)處 于熱平衡狀態(tài)為了計(jì)算簡(jiǎn)那么這個(gè)體積的形狀是無關(guān)緊要的。對(duì)于晶態(tài)的系統(tǒng)，元胞的形狀是有影響的。的情況下，則

12、其體積元胞的線度大小為L(zhǎng)便，對(duì)于氣體和液體，我們?nèi)∫粋€(gè)立方形的體積為MD元胞。設(shè)MD3。由于引進(jìn)這樣的立方體箱子，將產(chǎn)生六個(gè)我們不希望岀現(xiàn)的表面。模擬中碰撞這些箱的表面的 L為然而這些表面的存在對(duì)系統(tǒng)的任何一種粒子應(yīng)當(dāng)被反射回到元胞內(nèi)部，特別是對(duì)粒子數(shù)目很少的系統(tǒng)。我性質(zhì)都會(huì)有重大的影響。為了減小引入的表面效應(yīng)， 我們采用周期性邊界條件。采用這種邊界條件，這里我構(gòu)造岀一個(gè)準(zhǔn)無窮大的體積來更精確地代 表宏觀系統(tǒng)。實(shí)際上，們就可以消除引入的表面效應(yīng)，周期性邊界條件的數(shù)學(xué)表示們做了一個(gè)假定，即讓這個(gè)小體積元胞鑲嵌在一個(gè)無窮大的大塊物質(zhì)之中。形式為A（x） = A（x + »L）, 用也）（

13、3.5元胞完全等同地 MD、n。為任意整數(shù)。這個(gè)邊界條件就是命令基本其中A為任意的可觀測(cè)量，n、n3i2元胞的六方體表 MD重復(fù)無窮多次。具體在實(shí)現(xiàn)該邊界條件時(shí)是這樣操作的：當(dāng)有一個(gè) 粒子穿過基本 元胞內(nèi)。面時(shí)，就讓這個(gè)粒子以相同的速度穿過此表面對(duì)面的表面重新進(jìn)入該MD這個(gè)約定是在由無窮重復(fù)通常采用最小像力約定。在分子動(dòng)力學(xué)模擬中考慮粒子間的相互作用時(shí)，個(gè)粒子）中NN-1個(gè)（設(shè)在此元胞內(nèi)有的 MD基本元胞中，一個(gè)粒子只同它所在的基本元胞內(nèi)的 另外 之間的距離為處的粒子 jr處的粒子i同r的每個(gè)粒子或其最鄰近影像粒子發(fā)生相互作用。如果ji帀=皿冏彳+«£），(3.6)（對(duì)一切

14、的n ）的數(shù)值應(yīng)當(dāng)選得。通常LrcvL/2實(shí)際上這個(gè)約定就是通過滿足不等式條件來截?cái)辔粍?shì)（r為截止距離）c以避免有限尺寸效應(yīng)。采用最小像力約定使得在L/2的粒子的相互作用可以忽略，很大，使得距離大于-函數(shù)的奇異性，這會(huì)給模擬計(jì)算帶來誤差。截?cái)嗵幜Ｗ拥氖芰τ幸粋€(gè)S 2、給定初始條件不模擬進(jìn)人對(duì)系統(tǒng)微分方程組做數(shù)值求解的過程時(shí)，需要知道粒子的初始位置和速度的數(shù)值。MD方法需要兩組坐標(biāo)來啟動(dòng)計(jì)算；一組是零時(shí)刻的坐標(biāo)，另Verlet同的算法要求不同的初始條件。例如，一般來說系統(tǒng)的初始條件都是但是，一組是前進(jìn)一個(gè)時(shí)間步長(zhǎng)時(shí)的坐標(biāo)，或者是一組零時(shí)刻的速度值。因精確選擇待求系統(tǒng)的初始條件是沒有什么意義的，不

15、可能知道的。表面上看這是一個(gè)難題。實(shí)際上，常但是初始條件的合理選擇將可以加快系統(tǒng)趨于平衡。為模擬時(shí)間足夠長(zhǎng)時(shí)，系統(tǒng)就會(huì)忘掉初始條件。初始速度則從玻爾茲曼分布隨機(jī)令初始位置在差分劃分網(wǎng)格的格子上，1）用的初始條件可以選擇為： （）令初始位置隨機(jī)地偏 3）令初始位置隨機(jī)地偏離差分劃分網(wǎng)格 的格子，初始速度為零。（抽樣得到。（2離差分劃分網(wǎng)格的格子， 初始速度從玻爾茲曼分布隨機(jī) 抽樣得到。3 趨于平衡這樣計(jì)按照上面給岀的運(yùn)動(dòng)方程、邊界條件和初始條件， 就可以進(jìn)行分子動(dòng)力學(xué)模擬計(jì)算。但是，為了使系統(tǒng)達(dá)到平衡，算岀的系統(tǒng)不會(huì)具有所要求的系統(tǒng)能量，并且這個(gè)狀態(tài)本身也還不是一個(gè)平衡態(tài)。然直到系統(tǒng)具有所要求的

16、能量。模擬中需要一個(gè)趨衡過程。在這個(gè)過程中，增加或從系統(tǒng)中移出能量，我們稱；這時(shí)系統(tǒng)已經(jīng)達(dá)到使系統(tǒng)持續(xù)給出確定能量值。后，再對(duì)運(yùn)動(dòng)方程中的時(shí)間向前積分若于步，的大小選擇是十分重要 h平衡態(tài)。這段達(dá)到平衡所需的時(shí)間稱為弛豫時(shí)間。在MD模擬中，時(shí)間步長(zhǎng)必須取得小一些；但是取得太小，系統(tǒng)模擬的h的。它決定了模擬所需要的時(shí)間。為了減小誤差，步長(zhǎng)。例如，對(duì)一個(gè)具有幾百個(gè)氬h弛豫時(shí)間就很長(zhǎng)。這里需要積累一定的模擬經(jīng)驗(yàn)，選擇適當(dāng)?shù)臅r(shí)間步長(zhǎng)-2這里.量級(jí)，就可以得到好的相圖為L(zhǎng)ennard-Jones( Ar)分于的體系，如果采用位勢(shì)，我們發(fā)現(xiàn)取h10-14步，系統(tǒng)達(dá)10對(duì)應(yīng)的時(shí)間在1000秒的量級(jí)。如果模擬

17、hh選擇的是沒有量綱的，實(shí)際上這樣選擇的-11 10到平衡態(tài)，馳豫時(shí)間只有秒。44 宏觀物理量的計(jì)算實(shí)際計(jì)算宏觀物理量往往是在MD模擬的最后階段進(jìn)行的。它是沿著相空間軌跡求平均來計(jì)算得(N) A(0) r。如果已知初始位置和動(dòng)量為到的。例如對(duì)于一個(gè)宏觀物理量A，它的測(cè)量值應(yīng)當(dāng)為平均值(n)(0)(上標(biāo)N表示系綜Np個(gè)粒子的對(duì)應(yīng)坐標(biāo)和動(dòng)量參數(shù))，選擇某種MD算 法求解具有初值問題和(N)(N)(t)對(duì)軌跡平均的宏觀物理量 A的表示為(t) , p的運(yùn)動(dòng)方程，便得 到相空間軌跡rA - Iim(，/嚴(yán)的如P叫)(3.7)如果宏觀物理量為動(dòng)能，它的平均為lim下d叭("")0 )

18、 0(3.8)由于在模擬過程中計(jì)算出的動(dòng)能值是在不連續(xù)的路徑上的值，因此公式( 間的各個(gè)間斷點(diǎn)卩上計(jì)算動(dòng)能的平均值3.8 )可以表示為在時(shí)2m(3.9)在MD模擬過程中，溫度是需要加以監(jiān)測(cè)的物理量,特別是在模擬的起始階段。根據(jù)能量均分定理，我們可以從平均動(dòng)能值計(jì)算得到溫度值:(3.10)d= 3。系統(tǒng)內(nèi)部的位形能量的軌其中d為每個(gè)粒子的自由度，如果不考慮系統(tǒng)所受的約束，則歸池£»(護(hù))道平均值為(3.11)假定位勢(shì)在r處被截?cái)?，那么上式?jì)算岀的勢(shì)能以及由此得到的總能量就包含有誤差。為了對(duì)此偏差c作岀修正，我們采用對(duì)關(guān)聯(lián)函數(shù)來表示位能U J N 2矽 J i4(r)g(r)r

19、zdr(3.12)式中的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之間的平均粒子數(shù)，則g(JV n(r)=.N 4龍尸仏廠(3.13)在MD模擬過程中，所有的距離已經(jīng)在力的計(jì)算中得到，因而很容易計(jì)算對(duì)關(guān)聯(lián)函數(shù)的值。圖3.2為由計(jì)算機(jī)模擬得到的兩組不同參數(shù)下的對(duì)關(guān)聯(lián)函數(shù)的例子。由于位勢(shì)的截?cái)?，?duì)關(guān)聯(lián)函數(shù) 僅對(duì)rv L/2以c下的距離有意義。在公式(3.11 )中，所有的位能都加到截?cái)嗑嚯x為止，尾部修U c = 2兀p u(r)g(r)r

20、2dr正可以取為斥：(3.14)也可以利用含對(duì)壓強(qiáng)可以通過計(jì)算在面積元dA的法線方向上凈動(dòng)量轉(zhuǎn)移的時(shí)間平均值來得到,關(guān)聯(lián)函數(shù)的維里狀態(tài)方程計(jì)算。該維里狀態(tài)方程可以寫為5一項(xiàng)是對(duì)位一項(xiàng)是由相互作用力程之內(nèi)的貢獻(xiàn)引起的, 兩項(xiàng)， 勢(shì)截?cái)嗟母恼?xiàng)：g(r) 43r3dr)(3.15至于勢(shì)能的計(jì)算，我們可以把積分劃分為)(3.16g(廠)字4兀八"其中長(zhǎng)程改正項(xiàng)為：6 %dr)(3.17模擬：下面我們將討論具體如何進(jìn)行 MD0 L0 X0 3.0 4_0 5*00 r/or/<7由計(jì)算機(jī)模擬得到的兩組不同參數(shù)下的對(duì)關(guān)聯(lián)函數(shù)圖4、平衡態(tài)分子動(dòng)力學(xué)模擬模模擬。 兩種系統(tǒng)狀態(tài)的 MD 在經(jīng)典

21、 綜MD擬，另一類是對(duì)非平衡態(tài)的 NVT ( NVE )模擬，正則系綜的6.2模擬方法的應(yīng)用當(dāng)中, 微正則系綜的模擬。一種是對(duì)平衡態(tài)的MD模擬又可以分為如下類型：MD)NPHMD (模擬和等焓等壓系綜等溫等壓系綜MD ()模擬，MDMD存在著對(duì)對(duì)平衡態(tài)系面我們僅對(duì)平衡態(tài)的 MD MD、微正則系綜的模擬 模擬時(shí)，首先我們要確定所采用的相互作用模型。的多粒子體系，其粒子間的相互作用位勢(shì)是球?qū)ΨQ的,MD( NPT) 方法中前兩類模擬做簡(jiǎn)單的介紹。模擬等。下1)我們假定一個(gè)孤立在進(jìn)行對(duì)微正則系綜的 則其哈密頓量可以寫為MD個(gè)微正則系綜中，j 數(shù)也是不變的。此外； 系統(tǒng)受到的四個(gè)約束。)4.1(由于這個(gè)

22、系統(tǒng)的哈密頓量中不顯在這個(gè)粒子之間的距離。個(gè)粒子與第其中r是第ij由于整個(gè)系系統(tǒng)的體積和粒子式地岀現(xiàn)時(shí)間關(guān)聯(lián)，因而系統(tǒng)的能量是個(gè)守恒量。P恒等于零。這就是統(tǒng)并未運(yùn)動(dòng)，所以整個(gè)系統(tǒng)的總動(dòng)量由該系統(tǒng)的哈密頓量可以推導(dǎo)出牛頓方程形式的運(yùn)動(dòng)方程組= 2（勺）'（心 1,2 N）dr m '(4.2)的4.2Veriet要用數(shù)值求解的方法解岀(4.2)微分方程組，類似于本章第二節(jié)中介紹的方法方 程組(求解變成求解方程組：咽 + 防=2r(f) -譏+ Ft(f)h2 Im.。士人2N)(4.3)6該方程組反映岀：從前面 t和t-h時(shí)刻這兩步的空間坐標(biāo)位置及 t時(shí)刻的作用力，就可以算岀下

23、一步t+h時(shí)刻的坐標(biāo)位置。下面為了將(4.3)式寫成更簡(jiǎn)潔的形式，令屮)=巧山)，F(xiàn)# = Fg(4.4)則從(4.3)式可以得到如下差分方程組的形式r/fl+1) = 2r$)- 出宀 + Fh1 g (t =丄知(4.5)(o)(i)(2)(3) r, r r r,。，則通過求解方程組(4.5) 步步得到如果已知一組初始空間位置iiii由空間坐標(biāo)又可以算岀粒子的運(yùn)動(dòng)速度為(4.6)這里在第n+丨步算岀的速度是前一時(shí)刻；即第 n步的速度、因而動(dòng)能的計(jì)算比勢(shì)能的計(jì)算落 后一步。根據(jù)上述原理我們可以將粒子數(shù)恒定、體積恒定、能量恒定的微正則系綜(NVE )的MD模擬步驟設(shè)計(jì)如下：(0)(1) r,

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

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

26、 rhv2hm，計(jì)算在第n+ 1步時(shí)所有粒子所處的空間位 置（3）利用公式：iiii（n1）r 。i（n 1）（n 1）（n）（n 1）（n） vF/2（F vm h v 步時(shí)所有粒子的速度：n+1 （4）計(jì)算在第 訓(xùn)（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)能和位能；假如

27、總能量不等于給定恒定值，則通過對(duì)速度的調(diào)整來實(shí)現(xiàn)能量守恒。也就是將速度乘以一個(gè)標(biāo)度（seal ing）回子，該因子一般取為*11/2T*(N-D吩;4（4.10）然后再回到第一步，對(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模擬p Op） T和總動(dòng)量（在統(tǒng)計(jì)物理中的正則系綜模擬是針對(duì)一個(gè)粒子數(shù)N、體積V、溫度"為守恒量的系綜（NV

28、T ）。這種情況就如同一個(gè)系統(tǒng)置于熱浴之中，此時(shí)系統(tǒng)的能量可能有漲 落，但系統(tǒng)溫度則已經(jīng)保持恒定。在正則系綜的MD模擬中施加的約束與微正則系綜中的不一樣。正則系綜 MD方法是在運(yùn)動(dòng)方程組上加上動(dòng)能恒定（即溫度恒定）的約束，而不是像做正則系綜的MD模擬中對(duì)運(yùn)動(dòng)方程加上能量恒定的約束。在正則系綜MD的平衡化過程中，速度標(biāo)度因子一般選下面的形式較為合適（4.11）我們可將正則系綜 MD的Verlet算法的速度形式的模擬具體步驟列在下面:（1）t:（i1,2,N）?）給定初始空間位置（l i（1）v（2）給定初始速度i（n 1）5）（n）（n）2/hv 2r h rFm計(jì)算在第n+1 ）利用公式3（步

29、時(shí)所有粒子所處的空間位置iiii （n 1） r i 8）（nn1）n 1）（n）（2m） v F v h（F/，動(dòng)能和速度標(biāo)度因4（）計(jì)算在第n+ 1步時(shí)所有粒于的速度iiii子(3/V - 4)fcT1n v步粒子的速度值，并讓該值作為下一次計(jì)算時(shí)，第n+1 （5）計(jì)算將速度標(biāo)度因于Bi1 nn1 v v : ii步的模擬計(jì)算。n+ 2 （6）返回到步驟（3），開始第則退岀循環(huán)。這就是正則系綜的待系統(tǒng)達(dá)到平衡后，按照上面的步驟，對(duì)時(shí)間進(jìn)行一步步的循環(huán)模擬過程。MD模擬的應(yīng)用示例來看看模擬的結(jié)果。下面我們舉一個(gè)微正則系綜的MD模擬計(jì)算。對(duì)一個(gè)總能量確定的單原于（氬）粒子系統(tǒng)的MD 例：勢(shì)個(gè)原

30、子的模擬。粒于間的相互作用位勢(shì)為L(zhǎng)ennard -Jones我們具體選取2564.12 （）1/6處、該體系的粒于限制在一個(gè)立方體的箱br=2，其位置在&為能量單位）其中 -£為位勢(shì)的極小值（取1/22）/48和子，邊界上采用最小象力約定。我們采用自然單位制。長(zhǎng)度和時(shí)間的標(biāo) 度單位分別為b （m-12秒）10,這樣就使得運(yùn)動(dòng)方程為無（對(duì)氬原于該時(shí)間單位為3,它們分別具有）0.722 , 0.83134 p *）=（2.53 , 0.636）,（量綱形式。模擬時(shí)我們考慮兩個(gè)相圖上的點(diǎn)（T*，。初始條件假定為：各個(gè)原子處于一個(gè)面心立方格子的格點(diǎn)6.75L =兩種立方體的尺寸，即 L7.83和，用以比較其=3.6