【文章內(nèi)容簡介】 按照產(chǎn)生位形變化的方法，我們有兩類方法對有限的一系列態(tài)的物理量做統(tǒng)計平均：第一類是隨機模擬方法。它是實現(xiàn)Gibbs的統(tǒng)計力學途徑。在此方法中，體系位形的轉(zhuǎn)變是通過馬爾科夫(Markov)過程，由隨機性的演化引起的。這里的馬爾科夫過程相當于是內(nèi)稟動力學在概率方面的對應(yīng)物。該方法可以被用到?jīng)]有任何內(nèi)稟動力學模型體系的模擬上。隨機模擬方法計算的程序簡單，占內(nèi)存少，但是該方法難于處理非平衡態(tài)的問題。另一類為確定性模擬方法，即統(tǒng)計物理中的所謂分子動力學方法（Molecular Dynamics Method）。這種方法廣泛地用于研究經(jīng)典的多粒子體系的研究中。該方法是按該體系內(nèi)部的內(nèi)稟動力學規(guī)律來計算并確定位形的轉(zhuǎn)變。它首先需要建立一組分子的運動方程，并通過直接對系統(tǒng)中的一個個分子運動方程進行數(shù)值求解，得到每個時刻各個分子的坐標與動量，即在相空間的運動軌跡，再利用統(tǒng)計計算方法得到多體系統(tǒng)的靜態(tài)和動態(tài)特性, 從而得到系統(tǒng)的宏觀性質(zhì)。因此，分子動力學模擬方法可以看作是體系在一段時間內(nèi)的發(fā)展過程的模擬。在這樣的處理過程中我們可以看出：分子動力學方法中不存在任何隨機因素。系統(tǒng)的動力學機制決定運動方程的形式： 在分子動力學方法處理過程中，方程組的建立是通過對物理體系的微觀數(shù)學描述給出的。在這個微觀的物理體系中，每個分子都各自服從經(jīng)典的牛頓力學。每個分子運動的內(nèi)稟動力學是用理論力學上的哈密頓量或者拉格朗日量來描述，也可以直接用牛頓運動方程來描述。這種方法可以處理與時間有關(guān)的過程，因而可以處理非平衡態(tài)問題。但是使用該方法的程序較復(fù)雜，計算量大，占內(nèi)存也大。適用范圍廣泛： 原則上，分子動力學方法所適用的微觀物理體系并無什么限制。這個方法適用的體系既可以是少體系統(tǒng)，也可以是多體系統(tǒng)；既可以是點粒子體系，也可以是具有內(nèi)部結(jié)構(gòu)的體系；處理的微觀客體既可以是分子，也可以是其它的微觀粒子。 自五十年代中期開始，分子動力學方法得到了廣泛的應(yīng)用。它與蒙特卡洛方法一起已經(jīng)成為計算機模擬的重要方法。應(yīng)用分子動力學方法取得了許多重要成果，例如氣體或液體的狀態(tài)方程、相變問題、吸附問題等，以及非平衡過程的研究。其應(yīng)用已從化學反應(yīng)、生物學的蛋白質(zhì)到重離子碰撞等廣泛的學科研究領(lǐng)域。實際上，分子動力學模擬方法和隨機模擬方法一樣都面臨著兩個基本限制：一個是有限觀測時間的限制；另一個是有限系統(tǒng)大小的限制。通常人們感興趣的是體系在熱力學極限下（即粒子數(shù)目趨于無窮時）的性質(zhì)。但是計算機模擬允許的體系大小要比熱力學極限小得多，因此可能會出現(xiàn)有限尺寸效應(yīng)。為了減小有限尺寸效應(yīng)，人們往往引入周期性、全反射、漫反射等邊界條件。當然邊界條件的引入顯然會影響體系的某些性質(zhì)。數(shù)值求解時的離散化方法：對體系的分子運動方程組采用計算機進行數(shù)值求解時，需要將運動方程離散化為有限差分方程。常用的求解方法有歐拉法、龍格庫塔法、辛普生法等。數(shù)值計算的誤差階數(shù)顯然取決于所采用的數(shù)值求解方法的近似階數(shù)。原則上，只要計算機計算速度足夠大，內(nèi)存足夠多，我們可以使計算誤差足夠小。系綜的定義為：按照某些特定的量子態(tài)組成的集合體系，所謂的量子態(tài)就是P、T、V、N等狀態(tài)。它被用在分子動力學中是為了研究系統(tǒng)的微觀狀態(tài)與宏觀熱力學性質(zhì)的對應(yīng)規(guī)律。按照體系固定參數(shù)的不同，系綜可以分為多種，而每種系綜的配分函數(shù)也不盡相同。盡管在粒子數(shù)無限大的情況下，所有系綜是一樣的。這也就使宏觀性質(zhì)的計算式各不相同。下面就幾種常見的系綜及相應(yīng)的配風函數(shù)做簡單介紹。1，等溫等壓系綜一般情況下，實驗室所用的材料都認為是處在等溫等壓條件下的。等溫等壓系綜表示的是粒子數(shù)N、壓強P和體系溫度T固定的體系，簡寫為NPT。 其配分函數(shù)表示如下： ()式中，為普朗克常量為體系的哈密頓函數(shù)，和分別是坐標和動量。其Gibbs自由能表示為： ()式中，是波爾茲曼常數(shù)。NPT系綜是最常用的系綜，絕大多數(shù)的分子動力學模擬實驗都是在NPT系綜下進行的。在模擬研究中我們不僅要保證系綜溫度不變，還要保證系綜的壓強恒定，在模擬過程中，控制系綜的壓強是通過控制其體積來實現(xiàn)的。2，微正則系綜微正則系綜是指體系完全孤立，與外界不發(fā)生任何的物質(zhì)和能量交換，即粒子數(shù)N、體系總能量E和體系的總體積V都保持恒定，簡寫為NEV。配分函數(shù)表示如下： ()。微正則系綜的熱力學勢函數(shù)用熵來表示： ()在分子動力學模擬中，粒子的運動規(guī)律按經(jīng)典力學進行，因此，體系的能量是守恒的，在模擬體系的粒子數(shù)恒定，體積不變的情況下，模擬的系綜就是NEV系綜。3，正則系綜所謂正則系綜，就是體系的粒子總數(shù)N、總體積V和溫度T都保持恒定的體系，簡寫為NVT。其配分函數(shù)表示如下： () NVT系綜的熱力學勢函數(shù)為亥姆赫茲自由能，表示為： () 在分子動力學模擬中用到的系綜還有很多[51]，如具有確定粒子數(shù)N、壓強P和焓H的等壓等焓系綜（NPH），本文就不再一一介紹。 分子動力學基礎(chǔ)知識 基本理論分子動力學模擬的基本原理：分子動力學將連續(xù)介質(zhì)看成由N個原子或分子組成的粒子系統(tǒng)，各粒子之間的作用力可以通過量子力學勢能函數(shù)求導(dǎo)得出，忽略量子效應(yīng)后，運用經(jīng)典牛頓力學建立系統(tǒng)粒子運動數(shù)學模型,通過數(shù)值求解得到粒子在相空間的運動軌跡，然后由統(tǒng)計物理學原理得出該系統(tǒng)相應(yīng)的宏觀動態(tài)、靜態(tài)特性。分子動力學方法的出發(fā)點是對物理系統(tǒng)的確定的微觀描述。這個系統(tǒng)可以是一個少體系統(tǒng)，也可以是一個多體系統(tǒng)。這種描述可以是哈密頓描述或拉格朗日描述，也可以是直接用牛頓運動方程表示的描述。在前兩種情況下，運動方程必須應(yīng)用熟知的表述形式導(dǎo)出。顧名思義，分子動力學方法是用運動方程來計算系統(tǒng)的性質(zhì)，結(jié)果得到的既有系統(tǒng)的靜態(tài)特性，也有動態(tài)特性。對于由N個粒子構(gòu)成的系統(tǒng)，系統(tǒng)的總能量為動能與勢能之和。其中，系統(tǒng)動能由粒子的運動速度決定，假設(shè)模擬體系有n個原子構(gòu)成，用哈密頓運動方程表示的經(jīng)典力學的運動方程如下所示： () ()式中，i取，1,2…n，而哈密頓函數(shù)可表示為： ()式中，表示第i個粒子的廣義坐標，表示第i個粒子的動量，是該粒子的質(zhì)量，為系統(tǒng)的總勢能。根據(jù)牛頓運動方程： ()即可得到關(guān)于體系在各時刻的位置和速度。而其中的力是原子間相互作用勢來決定的，根據(jù)模擬計算中選用的勢函數(shù)，即可推導(dǎo)出體系中原子之間相互作用勢，即可計算出各個粒子在相空間中的運動軌跡，然后根據(jù)統(tǒng)計理論得到該系統(tǒng)的宏觀物理特性。分子動力學模擬是在原子尺度上模擬物質(zhì)的性質(zhì)，其根本問題是要確定一定量有相互作用的粒子的演化規(guī)律，即根據(jù)體系的動力學規(guī)律，計算出體系中每個粒子在每個時刻的坐標與速度，由此得到體系在相空間中的運動軌跡，然后運用統(tǒng)計學方法得到系統(tǒng)的靜態(tài)和動態(tài)特性，進而推倒系統(tǒng)的宏觀特性。本方法有一定的確定性，不存在任何的隨機過程，計算出的運動軌跡具有確切的動力學含義。分子動力學的基本假設(shè)主要有兩個，介紹如下：1，波恩奧本海默近似波恩奧本海默近似認為繞原子核運動的電子的運動速度比原子核的運動速度要大很多，因此可認為電子能夠及時的調(diào)整其自身位置，來適應(yīng)原子核位置的變化。因此，體系內(nèi)的相互作用就變成了一個只依賴兩個原子核間距離的勢函數(shù)的量，而電子的影響則被化為了一個平均的勢場。另外，原子核的質(zhì)量要比電子大很多，其運動規(guī)律可以用經(jīng)典力學理論，一般用牛頓方程來描述。2，各態(tài)歷經(jīng)假設(shè)各態(tài)歷經(jīng)假設(shè)認為，分子動力學由時間的平均來代替系綜平均以計算物質(zhì)的各種性能，即是用相空間中的一條軌跡而不是用整個相空間來確定物質(zhì)性質(zhì)，這種代替能行得通的前提是這條軌 能遍歷相空間中的所有狀態(tài)。 分子動力學計算中，必須要對一些方程求解，才能最終計算出體系中原子在各個時刻的坐標和速度。因此人們基于有限差分法，提出了一系列求解微分運動方程的算法，其中Verlet發(fā)展的算法在分子動力學中使用的最為普遍。 1967年，Verlet在研究LennardJones液體的性質(zhì)時提出了著名的Verlet算法，它是建立在體系中原子的位置、速度、加速度和Taylor展開式的基礎(chǔ)之上的： ()，可由t及的位置預(yù)測時的位置。進一步推導(dǎo)，可得到t時刻速度為： ()，t時的速度可由及的位置得到。，若很小，則容易導(dǎo)致誤差。 為了減小Verlet算法中的誤差，1970年，Hockney對其進行了修正，提出了新的算法，即所謂的蛙跳 (1eapflog)算法： （） （） （）雖然蛙跳算法修正了Verlet算法的誤差，使計算更準確，但是此算法計算出來的原子的位置和速度不同步，這就導(dǎo)致計算結(jié)果有一定的偏差，為了修正這種偏差，1982年，Swope等人[56]提出了速度Verlet，表示為： （） （） 速度Verlet算法解決了計算過程中原子位置和速度不同步的問題，更準確的表述了實際原子的狀態(tài)和性質(zhì)，而且Verlet算法的還有執(zhí)行簡明，對計算機內(nèi)存要求不高等優(yōu)點，因此被研究者廣泛使用。 分子動力學模擬的程序組成及邊界條件一．分子動力學模擬程序組成（即分子動力學模擬過程） 分子動力學方法的基本方程都是線性或非線性的二階常微分方程，對此人們很早就提出了阿達姆斯(Adams)方法、歐拉(Euler)折線法、改進的歐拉法及其龍格庫塔(RungeKutta)方法等，并在實際中得到了廣泛應(yīng)用。就分子動力學而言，主要使用Verlet法和Gear法進行數(shù)值積分求解。所謂分子動力學模擬，是指對于原子核和電子所構(gòu)成的多體系統(tǒng)，用計算機模擬原子核的運動過程，并從而計算系統(tǒng)的結(jié)構(gòu)和性質(zhì)，其中每一原子核被視為在全部其它原子核核電子所提供的經(jīng)驗勢場作用下按牛頓定律運動。分子動力學方法的出發(fā)點是對物理系統(tǒng)的確定的微觀描述。這個系統(tǒng)可以是一個少體系統(tǒng)，也可以是一個多體系統(tǒng)。這種描述可以是哈密頓描述或拉格朗日描述，也可以是直接用牛頓運動方程表示的描述。在前兩種情況下，運動方程必須應(yīng)用熟知的表述形式導(dǎo)出。顧名思義，分子動力學方法是用運動方程來計算系統(tǒng)的性質(zhì)，結(jié)果得到的既有系統(tǒng)的靜態(tài)特性，也有動態(tài)特性。而蒙特卡羅方法只能得到系統(tǒng)的位形特性，雖然它也有一個動力學解釋。 圖31 分子動力學的基本原理流程分子動力學模擬步驟如下：1，規(guī)定初始條件。令初始位置在格子的格點上，初始速度必須按照波爾茲曼分布得出。2，選擇適當?shù)膭莺瘮?shù)。所謂勢函數(shù)就是把體系中粒子之間相互作用表示出來。選擇適當?shù)膭莺瘮?shù)對模擬結(jié)果至關(guān)重要，本文采用的是原子嵌入勢(EAM)。3，弛豫平衡。按照上述方法建立的體系很可能不是處于平衡態(tài)，因此系統(tǒng)需要一個趨衡階段，即使系統(tǒng)弛豫若干時間步長，使系統(tǒng)達到平衡態(tài)。4，物理量的計算。按照系統(tǒng)輸入?yún)?shù)的設(shè)定，計算出相應(yīng)的計算結(jié)果。在計算機模擬出現(xiàn)以前，人們通常根據(jù)原子間相互作用力等微觀信息來推斷計算多原子體系的結(jié)構(gòu)和性質(zhì)。然而，原子相互作用力稍微復(fù)雜一些，不用說求解統(tǒng)計力學方程嚴格解，就是進行數(shù)值求解也是一件很困難的事情。分子動力學方法就是數(shù)值地分析求解多體系統(tǒng)的確定性運動方程，并根據(jù)對所求結(jié)果進行統(tǒng)計處理決定粒子的軌跡，從而給出物性預(yù)測和微觀結(jié)構(gòu)信息的一種模擬方法。隨著計算機的不斷發(fā)展，計算機軟件也得到了長足的進步，這使復(fù)雜體系的分子動力模擬計算成為可能。分子動力學模擬方法的計算機模擬程序也得到較大發(fā)展。一個完整的分子動力學模擬程序主要由以下幾部分組成：(1) 讀入指定運算條件的參數(shù)，如初始溫度、粒子數(shù)、密度和時間等。(2) 體系初始化，即選定初始位置和速度。(3) 計算作用于所有粒子上的力。這是幾乎所有的分子動力學模擬中最費時的部分。如果不使用技巧的話，力的計算時間與N2成正比。采用幾種有效的方法可以用來加速短程力和長程力的計算，這時計算時間與N成正比。(4) 求解牛頓運動方程。這一步和上一步構(gòu)成了模擬的核心。重復(fù)這兩步直至計算體系的演化到指定的時間長度。(5) 中心循環(huán)完成之后，計算并打印測定量的平均值，模擬結(jié)束。對于不同的系綜，粒子所遵循的運動規(guī)律不盡相同，因此，使用什么樣的分子運動方程需要具體問題具體分析。對于一個復(fù)雜的系綜，單個粒子所遵循的運動規(guī)律已經(jīng)完全確定，多種類型的粒子間相互作用的力的形式也完全確定，這樣對于一個系綜內(nèi)部粒子的運動的描述具有明確的表達公式，但是其運動方程的解析卻是不可能求得的。二．分子動力學的邊界條件對于不同的系綜，粒子所遵循的運動規(guī)律不盡相同，因此，使用什么樣的分子運動方程需要具體問題具體分析。對于一個復(fù)雜的系綜，單個粒子所遵循的運動規(guī)律已經(jīng)完全確定，多種類型的粒子間相互作用的力的形式也完全確定，這樣對于一個粒子系綜內(nèi)部粒子的運動的描述具有明確的表達公式，但是其運動方程的解析卻是不可能求得的。對于一小塊樣品，其內(nèi)部粒子數(shù)也是很多的，要想精確求出系綜中粒子間相互作用和運動狀態(tài)是不可能的，只有通過少量粒子的運動性質(zhì)計算來求出宏觀性質(zhì)。系綜總粒子數(shù)較少，相較于系綜總粒子，構(gòu)建的模擬計算晶胞的邊緣的粒子較多，邊緣粒子受力狀態(tài)與晶胞內(nèi)部粒子不同，為了降低