分子動(dòng)力學(xué)作業(yè)

1、分子動(dòng)力學(xué)(MD)1 分子動(dòng)力學(xué)（MD）基礎(chǔ) 1.1 MD分類(lèi) 1.2 MD簡(jiǎn)介 1.3 MD適用范圍2 分子動(dòng)力學(xué)運(yùn)動(dòng)方程數(shù)值求解 2.1 基礎(chǔ)知識(shí) 2.1.1 運(yùn)動(dòng)方程 2.1.2 空間描述 2.1.3 最小作用量原理 2.1.4 拉格朗日(Lagrange)方程 2.1.5 哈密頓(Hamilton)方程 2.2 粒子運(yùn)動(dòng)方程的數(shù)值解法 2.2.1 Verlet算法 2.2.2 歐拉(Euler)預(yù)測(cè)矯正公式 2.2.3 Gear預(yù)測(cè)矯正方法3 分子動(dòng)力學(xué)原胞與邊界條件 3.1 分子動(dòng)力學(xué)原胞 3.2 邊界條件 3.2.1 自由表面邊界 3.2.2 固定邊界 3.2.3 柔性邊界 3.2

2、.4 周期性邊界4 勢(shì)函數(shù)與分子力場(chǎng) 4.1 勢(shì)函數(shù) 4.1.1 兩體勢(shì) 4.1.2 多體勢(shì) 4.2 分子力場(chǎng) 4.2.1 分子力場(chǎng)函數(shù)的構(gòu)成 4.2.2 常用力場(chǎng)函數(shù)和分類(lèi)5 分子動(dòng)力學(xué)模擬的基本步驟 5.1 設(shè)定模擬所采用的模型 5.2 給定初始條件 5.3 趨于平衡計(jì)算 5.4 宏觀物理量的計(jì)算6 平衡態(tài)分子動(dòng)力學(xué)模擬 6.1 系綜 6.2 微正則系綜的分子動(dòng)力學(xué)模擬 6.3 正則系綜的分子動(dòng)力學(xué)模擬1 分子動(dòng)力學(xué)（MD）基礎(chǔ)1.1MD分類(lèi) 微正則系綜（VNE） 正則系綜（VNP） 平衡態(tài)MD 等溫等壓系綜（NPT）經(jīng)典MD 等焓等壓系綜（NPH） 巨正則系綜（） 非平衡態(tài)量子MD1.2

3、分子動(dòng)力學(xué)(MD)簡(jiǎn)介分子動(dòng)力學(xué)是在原子、分子水平上求解多體問(wèn)題的重要的計(jì)算機(jī)模擬方法。分子動(dòng)力學(xué)方法為確定性模擬方法，廣泛地用于研究經(jīng)典的多粒子體系的研究中,是按該體系內(nèi)部的內(nèi)稟動(dòng)力學(xué)規(guī)律來(lái)計(jì)算并確定位形的轉(zhuǎn)變。分子動(dòng)力學(xué)方法是通過(guò)建立一組分子的運(yùn)動(dòng)方程，并通過(guò)直接對(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ì)。在分子動(dòng)力學(xué)中，粒子的運(yùn)動(dòng)行為是通過(guò)經(jīng)典的Newton運(yùn)動(dòng)方程所描述。系統(tǒng)的所有粒子服從經(jīng)典力學(xué)的運(yùn)動(dòng)規(guī)律，它的動(dòng)力學(xué)方程就是從經(jīng)典力學(xué)的運(yùn)動(dòng)方程拉格朗日(la

4、grange)方程和哈密頓(Hamilton)方程導(dǎo)出。1.3適用范圍原則上，分子動(dòng)力學(xué)方法所適用的微觀物理體系并無(wú)什么限制。這個(gè)方法適用的體系既可以是少體系統(tǒng)，也可以是多體系統(tǒng)；既可以是點(diǎn)粒子體系，也可以是具有內(nèi)部結(jié)構(gòu)的體系；處理的微觀客體既可以是分子，也可以是其它的微觀粒子。實(shí)際上，分子動(dòng)力學(xué)模擬方法和隨機(jī)模擬方法一樣都面臨著兩個(gè)基本限制：一個(gè)是有限觀測(cè)時(shí)間的限制；另一個(gè)是有限系統(tǒng)大小的限制。通常人們感興趣的是體系在熱力學(xué)極限下（即粒子數(shù)目趨于無(wú)窮時(shí)）的性質(zhì)。但是計(jì)算機(jī)模擬允許的體系大小要比熱力學(xué)極限小得多，因此可能會(huì)出現(xiàn)有限尺寸效應(yīng)。為了減小有限尺寸效應(yīng)，人們往往引入周期性、全反射、漫反

5、射等邊界條件。當(dāng)然邊界條件的引入顯然會(huì)影響體系的某些性質(zhì)。2 分子動(dòng)力學(xué)運(yùn)動(dòng)方程數(shù)值求解2.1 基礎(chǔ)知識(shí)2.1.1 運(yùn)動(dòng)方程系統(tǒng)的動(dòng)力學(xué)機(jī)制決定運(yùn)動(dòng)方程的形式。在分子動(dòng)力學(xué)方法處理過(guò)程中，方程組的建立是通過(guò)對(duì)物理體系的微觀數(shù)學(xué)描述給出的。在這個(gè)微觀的物理體系中，每個(gè)分子都各自服從經(jīng)典的牛頓力學(xué)。每個(gè)分子運(yùn)動(dòng)的內(nèi)稟動(dòng)力學(xué)是用理論力學(xué)上的哈密頓量或者拉格朗日量來(lái)描述，也可以直接用牛頓運(yùn)動(dòng)方程來(lái)描述。采用分子動(dòng)力學(xué)方法時(shí)，必須對(duì)一組分子運(yùn)動(dòng)微分方程做數(shù)值求解。從計(jì)算數(shù)學(xué)的角度來(lái)看，這是個(gè)求一個(gè)初值問(wèn)題的微分方程的解。實(shí)際上計(jì)算數(shù)學(xué)為了求解這種問(wèn)題已經(jīng)發(fā)展了許多的算法。 但是并不是所有的這些算法都可以

6、用來(lái)解決物理問(wèn)題。2.1.2 空間描述在空間描述如何物體的運(yùn)動(dòng)，如果其本身的大小可以忽略時(shí)，就可以將其看作是粒子（或質(zhì)點(diǎn)）。粒子描述：空間位置：r 速度：v = dr/dt 加速度： 若一個(gè)系統(tǒng)由N個(gè)粒子組成，則粒子描述： 空間位置：r1,r2,r3,rN 笛卡爾坐標(biāo)系，粒子有3N個(gè)自由度設(shè)系統(tǒng)有s個(gè)自由度 廣義坐標(biāo)：q1,q2,q3,qN 廣義速度：q1,q2,q3,qN2.1.3 最小作用量原理莫培督1744年提出最小作用量原理：保守的、完整的力學(xué)系統(tǒng)，由某一初位形轉(zhuǎn)變到另一位形的一切具有相同能量的可能運(yùn)動(dòng)中，真實(shí)的運(yùn)動(dòng)是其作用量具有極小值的那種運(yùn)動(dòng)。力學(xué)系統(tǒng)中，構(gòu)造能量函數(shù)L及其作用量S

7、作用量的積分式叫做泛函(functional)，作用量取極值的方法就是求其變分S = 0。2.1.4 拉格朗日(Lagrange)方程由最小作用量原理可導(dǎo)出拉格朗日方程對(duì)于孤立的保守系統(tǒng)，每個(gè)粒子在勢(shì)場(chǎng)U中運(yùn)動(dòng)，則系統(tǒng)整體的Lagrange函數(shù)是得到第i個(gè)粒子的牛頓運(yùn)動(dòng)方程（指每個(gè)粒子的自由度）2.1.5 哈密頓(Hamilton)方程哈密頓(Hamilton)原理:保守的、完整的力學(xué)系統(tǒng)在相同時(shí)間內(nèi)，由某一初位形轉(zhuǎn)移到另一已知位形的一切可能運(yùn)動(dòng)中，真實(shí)運(yùn)動(dòng)的作用函數(shù)具有極值，即作用函數(shù)的變分等于零。哈密頓(Hamilton)方程，Lagrange函數(shù)全微分形式： 則定義哈密頓函數(shù)或哈密頓量為

8、：哈密頓函數(shù)H是動(dòng)量和坐標(biāo)的函數(shù)，是動(dòng)能和勢(shì)能之和：變量為動(dòng)量p和坐標(biāo)r的Hamilton方程：這就是變量為動(dòng)量p和坐標(biāo)q的哈密頓方程。如果系統(tǒng)的哈密頓函數(shù)不顯含時(shí)間，就有dH/dt=0，即得到能量守恒定律。 2.2 粒子運(yùn)動(dòng)方程的數(shù)值解法設(shè)粒子的坐標(biāo)、速度、動(dòng)量及其作用力分別用x(t)，v(t)，p(t)，f(x,t)表示，其初始值為x(0)，v(0)，p(0)，f(0)。則決定粒子運(yùn)動(dòng)的牛頓方程是同時(shí)由于數(shù)值計(jì)算求解微分方程是用差分的方法，習(xí)慣上將時(shí)間的變化間隔 t 用 h 表示，叫做時(shí)間步長(zhǎng)。坐標(biāo)預(yù)測(cè)的Taylor展開(kāi)式： 勢(shì)函數(shù) U 應(yīng)看作是離子間相互作用勢(shì) V 和外勢(shì) Uext2.2

9、.1 Verlet算法根據(jù)Taylor公式，在t 時(shí)刻求t +h 時(shí)刻的坐標(biāo)和作用力時(shí)，分成向前和向后的Taylor展開(kāi)式將上面兩式相加得到從而得到坐標(biāo)的計(jì)算公式將兩個(gè)Taylor展開(kāi)式相減可以得到從而得到速度或動(dòng)量的計(jì)算公式和加速度或作用力的計(jì)算公式于是得到用上一個(gè)時(shí)間點(diǎn)(t-h)和當(dāng)前時(shí)間點(diǎn)(t)的坐標(biāo)r和速度v及作用力f，來(lái)計(jì)算出下一個(gè)時(shí)間點(diǎn)(t+h)的坐標(biāo)和速度及作用力的三點(diǎn)公式。這種方法稱(chēng)為Verlet方法。2.2.2 歐拉(Euler)預(yù)測(cè)矯正公式第一種方法，開(kāi)放式或單預(yù)測(cè)式。對(duì)下一步的位置和力的計(jì)算僅與前幾步的已知位置和力有關(guān)。實(shí)際計(jì)算證明，這種計(jì)算方式的精度小而誤差大，所以實(shí)際

10、計(jì)算中除了在特殊情況下，已不再采用。第二種方法，封閉式或預(yù)測(cè)矯正式。它的計(jì)算步驟是：用預(yù)測(cè)公式計(jì)算得到下一步的數(shù)值yn+1；將其帶入函數(shù)式，計(jì)算f(yn+1)去矯正預(yù)測(cè)值yn+1；從而得到矯正后的rn+1。具體操作看下面的歐拉(Euler)預(yù)測(cè)矯正公式：預(yù)測(cè)值矯正值2.2.3 Gear預(yù)測(cè)矯正方法Gear發(fā)展出預(yù)測(cè)矯正方法(Predict-corrector)。經(jīng)證明，這是一種精度很高的完全適用于分子動(dòng)力學(xué)的算法被廣泛應(yīng)用。為方便，使用矢量記法。將下一步預(yù)測(cè)值的每一項(xiàng)進(jìn)行Taylor展開(kāi)組成一個(gè)列矢量，為書(shū)寫(xiě)方便記做轉(zhuǎn)置形式稱(chēng)為N-表象矢量。若將當(dāng)前和以前幾個(gè)時(shí)刻的坐標(biāo)值作為一個(gè)矢量的各元素，

11、則有稱(chēng)為C-表象矢量。還有一種使用起來(lái)比較方便的F-表象，矢量的元素由當(dāng)前的坐標(biāo)、速度和當(dāng)前幾前兩步的力（或加速度）構(gòu)成稱(chēng)為F-表象矢量。Gear指出，各種表象可以通過(guò)一個(gè)變換進(jìn)行相互轉(zhuǎn)換。矢量rn+1也就是所有n-1階多項(xiàng)式空間中的矢量，將一個(gè)表象R1變換到另一個(gè)表象R2的變換完全是一個(gè)在該空間中的正交線(xiàn)形變換其中，變換矩陣T可以很容易找到。(1) 預(yù)測(cè)算法。用矩陣代數(shù)的形式描述預(yù)測(cè)值為：yn+1是rn+1的預(yù)測(cè)值。 預(yù)測(cè)矩陣A是從某一物理量的對(duì)時(shí)間步長(zhǎng)的Taylor展開(kāi)中得到。N-表象中，5值預(yù)測(cè)矩陣為： A中的各列元素值構(gòu)成從0階到5階的二項(xiàng)式的展開(kāi)系數(shù)。更高階的預(yù)測(cè)矩陣只需在矩陣最后一

12、列添加更高階二項(xiàng)式的展開(kāi)系數(shù)。預(yù)測(cè)矩陣的變換矩陣為由此不難得到預(yù)測(cè)矩陣在其他表象中的形式。一個(gè)4值F-表象中的預(yù)測(cè)矩陣為(2) 矯正算法。 為了得到r 的下一步預(yù)測(cè)值rk+1，需要將矢量rk+1的每一項(xiàng)進(jìn)行Taylor展開(kāi)，由于t=kh，得到rk+1=B rk，B是(n+1)×(n+1)的系數(shù)矩陣，其第i 行第j 列的元素值為把預(yù)測(cè)項(xiàng)和矯正項(xiàng)結(jié)合起來(lái)C是矯正系數(shù)矢量。通過(guò)求解上面方程的本征值，即可求出系數(shù)向量C的各元素值。如果Taylor展開(kāi)的項(xiàng)數(shù)不太多，C的各值也可以直接用代入法求得。Gear推導(dǎo)的48值N-表象中矯正矢量C的值為代人預(yù)測(cè)值yn+1，計(jì)算預(yù)測(cè)位置上的加速度f(wàn)(yn+

13、1)/m，得到更為準(zhǔn)確的下一步位置可以由矯正式：3 分子動(dòng)力學(xué)原胞與邊界條件3.1 分子動(dòng)力學(xué)原胞分子動(dòng)力學(xué)模擬方法往往用于研究大塊物質(zhì)在給定密度下的性質(zhì)，而實(shí)際計(jì)算模擬不可能在幾乎是無(wú)窮大的系統(tǒng)中進(jìn)行。所以必須引進(jìn)一個(gè)叫做分子動(dòng)力學(xué)原胞的體積元，以維持一個(gè)恒定的密度。對(duì)氣體和液體，如果所占體積足夠大，并且系統(tǒng)處于熱平衡狀態(tài)的情況下，那么這個(gè)體積的形狀是無(wú)關(guān)緊要的。對(duì)于晶態(tài)的系統(tǒng)，原胞的形狀是有影響的。為了計(jì)算簡(jiǎn)便，取一個(gè)立方形的體積為分子動(dòng)力學(xué)原胞。設(shè)原胞的線(xiàn)度大小為L(zhǎng)，則體積為L(zhǎng)3。由于引進(jìn)這樣的立方體箱子，將產(chǎn)生六個(gè)我們不希望出現(xiàn)的表面。模擬中碰撞這些箱子的表面的粒子應(yīng)當(dāng)被反射回到原胞內(nèi)

14、部，特別是對(duì)粒子數(shù)目很少的系統(tǒng)。然而這些表面的存在對(duì)系統(tǒng)的任何一種性質(zhì)都會(huì)有重大的影響。3.2 邊界條件采用分子動(dòng)力學(xué)方法，必需對(duì)被計(jì)算的粒子系統(tǒng)給定適當(dāng)?shù)倪吔鐥l件。這些邊界條件大致可分成四種。（1） 自由表面邊界(free-surface boundary)這種邊界條件常用于大型的自由分子的模擬。（2） 固定邊界(rigid boundary)在所有要計(jì)算到粒子晶胞之外還要包上幾層結(jié)構(gòu)相同的位置不變的粒子，包層的厚度必須大于粒子間相互作用的力程范圍。包層部分代表了與運(yùn)動(dòng)粒子起作用的宏觀晶體的那一部分。 （3） 柔性邊界 (flexible boundary)這種邊界比固定邊界更接近實(shí)際。它允

15、許邊界上的粒子有微小的移動(dòng)以反映內(nèi)層粒子的作用力施加到它們身上時(shí)的情況。 （4） 周期性邊界(periodic boundary)在模擬較大的系統(tǒng)時(shí)，為了消除表面效應(yīng)或邊界效應(yīng)，常采用周期性邊界條件。就是讓原胞上下、左右、前后對(duì)邊上的粒子間有相互作用。 為了將分子動(dòng)力學(xué)原胞有限立方體內(nèi)的模擬，擴(kuò)展到真實(shí)大系統(tǒng)的模擬，通常采用周期性邊界條件。采用周期性邊界條件，就可以消除引入原胞后的表面效應(yīng)，構(gòu)造出一個(gè)準(zhǔn)無(wú)窮大的體積來(lái)更精確地代表宏觀系統(tǒng)。實(shí)際上，這里我們做了一個(gè)假定，即讓這個(gè)小體積原胞鑲嵌在一個(gè)無(wú)窮大的大塊物質(zhì)之中。做一個(gè)假設(shè)：讓這個(gè)小體積原胞鑲嵌在一個(gè)無(wú)窮大的大塊物質(zhì)中，周期性邊界條件的數(shù)學(xué)

16、表示形式為：其中A為任意可觀測(cè)量，n1,n2,n3為任意整數(shù)，這個(gè)邊界條件就是命令基本MD原胞完全等同地重復(fù)無(wú)窮多次。4 勢(shì)函數(shù)與分子力場(chǎng)4.1 勢(shì)函數(shù)由于物質(zhì)系統(tǒng)的復(fù)雜性以及原子間相互作用類(lèi)型的不同，很難得到滿(mǎn)足各種不同體系和物質(zhì)的一般性而又精度很高的勢(shì)函數(shù)，所以針對(duì)不同的物質(zhì)體系人們陸續(xù)發(fā)展了大量的經(jīng)驗(yàn)和半經(jīng)驗(yàn)的勢(shì)函數(shù)。半經(jīng)驗(yàn)勢(shì)：形式經(jīng)過(guò)一定的理論分析而得到，但其中一些參數(shù)則需要根據(jù)宏觀實(shí)驗(yàn)參數(shù)用經(jīng)驗(yàn)方法來(lái)確定，這些宏觀實(shí)驗(yàn)參數(shù)主要有彈性常數(shù)、平衡點(diǎn)陣常數(shù)以及內(nèi)聚能、空位形成能和層錯(cuò)能等。經(jīng)驗(yàn)勢(shì)：為了擬合的方便，在選擇勢(shì)函數(shù)的形式時(shí)，并不一定要求有確切的理論依據(jù)，而是出于經(jīng)驗(yàn)的估計(jì)和擬合方

17、便的需要，相對(duì)自由的選擇勢(shì)函數(shù)的形式。勢(shì)函數(shù)一般分為二體勢(shì)能和三體以上的多體勢(shì)能函數(shù)。勢(shì)函數(shù)形式：4.1.1 兩體勢(shì)兩體勢(shì)認(rèn)為原子之間的相互作用是兩兩原子之間的作用，與其他原子的位置無(wú)關(guān)。這種勢(shì)的形式可以比較充分地描述除半導(dǎo)體、金屬以外的無(wú)機(jī)化合物材料中的相互作用。（1）硬球模型（2）Lennard-Jones(LJ)勢(shì) （3）Morse勢(shì) （4）Borna-Mayer勢(shì) 4.1.2 多體勢(shì)兩體勢(shì)模型在解決分子晶體和離子型化合物問(wèn)題中取得了成功，但對(duì)于含有共價(jià)鍵的過(guò)渡金屬卻遇到了很多困難。（1）嵌入原子勢(shì)（EAM，Embedded atom model）（2）修正的嵌入原子勢(shì)（MEAM勢(shì)） 以

18、上討論的勢(shì)都是中心對(duì)稱(chēng)的，對(duì)于由共價(jià)鍵結(jié)合的有機(jī)分子和半導(dǎo)體材料是不適宜的。為了更好的描述各種含有共價(jià)鍵作用的材料，發(fā)展了很多考慮角度效應(yīng)的多體勢(shì)，其中修正的嵌入原子勢(shì)就是比較常用的一種，它的形式與嵌入原子勢(shì)一樣，只是在背景電子密度中考慮了角度效應(yīng)，其形式為： （3）Stillinger-Weber勢(shì) 4.2分子力場(chǎng) 根據(jù)量子力學(xué)的波恩-奧本海默近似，一個(gè)分子的能量可以近似看作構(gòu)成分子的各個(gè)原子的空間坐標(biāo)的函數(shù)。簡(jiǎn)單地講就是分子的能量隨分子構(gòu)型的變化而變化，描述這種分子能量和分子結(jié)構(gòu)之間關(guān)系的就是分子力場(chǎng)函數(shù)。分子力場(chǎng)函數(shù)為來(lái)自實(shí)驗(yàn)結(jié)果的經(jīng)驗(yàn)公式，可以講對(duì)分子能量的模擬比較粗糙，但是相比于精確

19、的量子力學(xué)從頭計(jì)算方法，分子力場(chǎng)方法的計(jì)算量要小數(shù)十倍。而且在適當(dāng)?shù)姆秶鷥?nèi)，分子力場(chǎng)方法的計(jì)算精度與量子化學(xué)計(jì)算相差無(wú)幾，因此對(duì)大分子復(fù)雜體系而言，分子力場(chǎng)方法是一套行之有效的方法。以分子力場(chǎng)為基礎(chǔ)的分子力學(xué)計(jì)算方法在分子動(dòng)力學(xué)、蒙特卡羅、分子對(duì)接等分子模擬方法中有著廣泛的應(yīng)用。4.2.1 分子力場(chǎng)函數(shù)的構(gòu)成 一般而言，分子力場(chǎng)函數(shù)由以下幾個(gè)部分構(gòu)成：鍵伸縮能：構(gòu)成分子的各個(gè)化學(xué)鍵在鍵軸方向上的伸縮運(yùn)動(dòng)所引起的能量變化。鍵角彎曲能：鍵角變化引起的分子能量變化。二面角扭曲能：?jiǎn)捂I旋轉(zhuǎn)引起分子骨架扭曲所產(chǎn)生的能量變化。非鍵相互作用：包括范德華力、靜電相互作用等與能量有關(guān)的非鍵相互作用。交叉能量項(xiàng)：

20、上述作用之間耦合引起的能量變化。構(gòu)成一套力場(chǎng)函數(shù)體系需要有一套聯(lián)系分子能量和構(gòu)型的函數(shù)；還需要給出各種不同原子在不同成鍵狀況下的物理參數(shù)，比如正常的鍵長(zhǎng)、鍵角、二面角等；這些力場(chǎng)參數(shù)多來(lái)自實(shí)驗(yàn)或者量子化學(xué)計(jì)算。勢(shì)函數(shù)形式, 但一般總表達(dá)為分子內(nèi)與分子間勢(shì)能之和:分子內(nèi)勢(shì)能(鍵合)包括鍵伸縮、鍵角彎曲和二面角扭轉(zhuǎn)勢(shì)能；分子間勢(shì)能(非鍵合)包括范德華勢(shì)和靜電勢(shì), 有的還包括H鍵：鍵合勢(shì)函數(shù)中，一些力場(chǎng)還包含交叉項(xiàng)，使精度更高。交叉項(xiàng)的含義：如鍵長(zhǎng)變化時(shí)，鍵角彎曲勢(shì)能隨鍵長(zhǎng)的不同而不同。4.2.2 常用力場(chǎng)函數(shù)和分類(lèi)不同的分子力場(chǎng)會(huì)選取不同的函數(shù)形式來(lái)描述能量與體系構(gòu)型之間的關(guān)系。到目前，不同的科研

21、團(tuán)隊(duì)設(shè)計(jì)了很多適用于不同體系的力場(chǎng)函數(shù)，根據(jù)他們選擇的函數(shù)和力場(chǎng)參數(shù)，可以分為以下幾類(lèi)：傳統(tǒng)力場(chǎng)第二代力場(chǎng)-第二代的勢(shì)能函數(shù)形式比傳統(tǒng)力場(chǎng)要更加復(fù)雜，涉及的力場(chǎng)參數(shù)更多，計(jì)算量也更大，當(dāng)然也相應(yīng)地更加準(zhǔn)確。通用力場(chǎng)-通用力場(chǎng)也叫基于規(guī)則的力場(chǎng)，它所應(yīng)用的力場(chǎng)參數(shù)是基于原子性質(zhì)計(jì)算所得，用戶(hù)可以通過(guò)自主設(shè)定一系列分子作為訓(xùn)練集來(lái)生成合用的力場(chǎng)參數(shù)。5 分子動(dòng)力學(xué)模擬的基本步驟5.1 設(shè)定模擬所采用的模型模型的設(shè)定，也就是勢(shì)函數(shù)的選取。勢(shì)函數(shù)的研究和物理系統(tǒng)上對(duì)物質(zhì)的描述研究息息相關(guān)。設(shè)定模型是分子動(dòng)力學(xué)模擬的第一步工作。例如在一個(gè)分子系統(tǒng)中，假定兩個(gè)分子間的相互作用勢(shì)為硬球勢(shì)，其勢(shì)函數(shù)表示為常用

22、的是Lennard-Jones 勢(shì)函數(shù)、EAM勢(shì)函數(shù)等，不同的物質(zhì)狀態(tài)描述用不同的勢(shì)函數(shù)。5.2 給定初始條件運(yùn)動(dòng)方程的求解需要知道粒子的初始位置和速度，不同的算法要求不同的初始條件。一般講，系統(tǒng)的初始條件不可能知道，實(shí)際上也不需要精確選擇代求系統(tǒng)的初始條件，因?yàn)槟M實(shí)踐足夠長(zhǎng)時(shí)，系統(tǒng)就會(huì)忘掉初始條件。當(dāng)然，合理的初始條件可以加快系統(tǒng)趨于平衡的時(shí)間和步伐，獲得好的精度。常用的初始條件可以選擇為：令初始位置在差分劃分網(wǎng)格的格子上，初始速度則從玻爾茲曼分布隨機(jī)抽樣得到；令初始位置隨機(jī)的偏離差分劃分網(wǎng)格的格子上，初始速度為零；令初始位置隨機(jī)的偏離差分劃分網(wǎng)格的格子上，初始速度也是從玻爾茲曼分布隨機(jī)抽

23、樣得到。5.3 趨于平衡計(jì)算按照給出的運(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è)趨衡過(guò)程。在這個(gè)過(guò)程中，增加或從系統(tǒng)中移出能量，直到系統(tǒng)具有所要求的能量。然后，再對(duì)運(yùn)動(dòng)方程中的時(shí)間向前積分若干步，使系統(tǒng)持續(xù)給出確定能量值。稱(chēng)這時(shí)系統(tǒng)已經(jīng)達(dá)到平衡態(tài)。這段達(dá)到平衡所需的時(shí)間稱(chēng)為弛豫時(shí)間。分子動(dòng)力學(xué)計(jì)算的基本思想是：賦予分子體系初始運(yùn)動(dòng)狀態(tài)之后，利用分子的自然運(yùn)動(dòng)在相空間中抽取樣本進(jìn)行統(tǒng)計(jì)分子動(dòng)力學(xué)計(jì)算，時(shí)間步長(zhǎng)就是抽樣的間隔，因而時(shí)間步長(zhǎng)的選取對(duì)動(dòng)力學(xué)模擬非常重要

24、。它決定了模擬所需要的時(shí)間。時(shí)間步長(zhǎng) h太長(zhǎng)會(huì)造成分子間激烈碰撞，體系數(shù)據(jù)溢出；為了減小誤差，步長(zhǎng) h 必須取得小一些；但是取得太小，系統(tǒng)模擬的弛豫時(shí)間就很長(zhǎng)，太短的時(shí)間步長(zhǎng)還會(huì)降低模擬過(guò)程搜索相空間的能力。因此一般選取的時(shí)間步長(zhǎng)為體系各個(gè)自由度中最短運(yùn)動(dòng)周期的十分之一。這里需要積累一定的模擬經(jīng)驗(yàn)，選擇適當(dāng)?shù)臅r(shí)間步長(zhǎng)h 。5.4 宏觀物理量的計(jì)算實(shí)際計(jì)算宏觀物理量往往是在分子動(dòng)力學(xué)模擬的最后階段進(jìn)行的。它是沿著相空間軌跡求平均來(lái)計(jì)算得到的。如，在模擬過(guò)程中計(jì)算出的動(dòng)能值是在不連續(xù)的路徑上的值，在時(shí)間的各個(gè)間斷點(diǎn)上計(jì)算動(dòng)能的平均值在分子動(dòng)力學(xué)模擬過(guò)程中，溫度是需要加以監(jiān)測(cè)的物理量，特別是在模擬的

25、起始階段。根據(jù)能量均分定理，我們可以從平均動(dòng)能值計(jì)算得到溫度值，其中d為每個(gè)粒子的自由度，如果不考慮系統(tǒng)所受的約束，則d3。6 平衡態(tài)分子動(dòng)力學(xué)模擬在經(jīng)典分子動(dòng)力學(xué)模擬方法的應(yīng)用當(dāng)中，存在著對(duì)兩種系統(tǒng)狀態(tài)的分子動(dòng)力學(xué)模擬。一種是對(duì)平衡態(tài)的MD模擬；另一類(lèi)是對(duì)非平衡態(tài)的分子動(dòng)力學(xué)模擬。 對(duì)平衡態(tài)系綜分子動(dòng)力學(xué)模擬又可以分為如下類(lèi)型： 微正則系綜的分子動(dòng)力學(xué)(N,V,E）模擬 正則系綜的分子動(dòng)力學(xué)（N,V,T）模擬 等溫等壓系綜分子動(dòng)力學(xué)(N,P,T) 模擬 等焓等壓系綜分子動(dòng)力學(xué)(N,P,H) 模擬巨正則系綜分子動(dòng)力學(xué)(V,T, ) 模擬6.1 系綜系綜(ensemble) 是指在一定的宏觀條件

26、下(約束條件)，大量性質(zhì)和結(jié)構(gòu)完全相同的、處于各種運(yùn)動(dòng)狀態(tài)的、各自獨(dú)立的系統(tǒng)的集合。全稱(chēng)為統(tǒng)計(jì)系綜。系綜是用統(tǒng)計(jì)方法描述熱力學(xué)系統(tǒng)的統(tǒng)計(jì)規(guī)律性時(shí)引入的一個(gè)基本概念；系綜是統(tǒng)計(jì)理論的一種表述方式，系綜理論使統(tǒng)計(jì)物理成為普遍的微觀統(tǒng)計(jì)理論 ；系綜并不是實(shí)際的物體，構(gòu)成系綜的系統(tǒng)才是實(shí)際物體。約束條件是由一組外加宏觀參量來(lái)表示。在平衡統(tǒng)計(jì)力學(xué)范疇下，可以用來(lái)處理穩(wěn)定系綜。根據(jù)宏觀約束條件，系綜被分為以下幾種：（1）正則系綜(canonical ensemble)，全稱(chēng)應(yīng)為“宏觀正則系綜”，簡(jiǎn)寫(xiě)為NVT，即表示具有確定的粒子數(shù)(N)、體積(V)、溫度(T)。正則系綜是蒙特卡羅方法模擬處理的典型代表。假

27、定N個(gè)粒子處在體積為V的盒子內(nèi)，將其埋入溫度恒為T(mén)的熱浴中。此時(shí)，總能量(E)和系統(tǒng)壓強(qiáng)(P)可能在某一平均值附近起伏變化。平衡體系代表封閉系統(tǒng)，與大熱源熱接觸平衡的恒溫系統(tǒng)。正則系綜的特征函數(shù)是亥姆霍茲自由能F(N,V,T)。（2） 微正則系綜(micro-canonical ensemble)，簡(jiǎn)寫(xiě)為NVE，即表示具有確定的粒子數(shù)(N)、體積(V)、總能量(E)。微正則系綜廣泛被應(yīng)用在分子動(dòng)力學(xué)模擬中。假定N個(gè)粒子處在體積為V的盒子內(nèi)，并固定總能量(E)。此時(shí)，系綜的溫度(T)和系統(tǒng)壓強(qiáng)(P)可能在某一平均值附近起伏變化。平衡體系為孤立系統(tǒng)，與外界即無(wú)能量交換，也無(wú)粒子交換。微正則系綜的特征函數(shù)是S(N,V,E)。（3）等溫等壓系綜(constant-pressure,constant-temperature)，簡(jiǎn)寫(xiě)為NPT，即表示具有確定的粒子數(shù)(N)、壓強(qiáng)(P)