SILICONLLEY■■【高新技術(shù)產(chǎn)業(yè)發(fā)展分子動力學(xué)模擬方法概述周曉平田壯壯忽曉偉(鄭州大學(xué)西亞斯國際學(xué)院河南鄭州451100)摘要:主要介紹分子動力學(xué)模擬的基本原理,闡述分子動力學(xué)方法的運動方程、數(shù)值解法、勢函數(shù)、邊界條件、適用系綜以及體系相關(guān)性質(zhì)的計算。最后指出分子動力學(xué)模擬方法的優(yōu)勢和發(fā)展方向。關(guān)鍵詞:分子動力學(xué);勢函數(shù);邊界條件中圖分類號:0414文獻標(biāo)識碼:A文章編號:1671-7597(2012)1210040-021分子動力模擬計算的基本原理利用蛙跳法計算僅需儲存(-)與r()兩種資料,可節(jié)分子動力計算的基本原理,即為利用牛頓運動定律。在分子動力模擬中,體系原子的一系列位移是通過對牛頓運動方程省儲存空間。其缺點是位置與速度不同步。這意味著在位置的積分得到的,結(jié)果是一條運動軌跡,它表明了系統(tǒng)內(nèi)原子的定時,不可能同時計算動能對總動能的貢獻位置與速度是如何隨時間而發(fā)生變化2.2Gea算法[1先由系統(tǒng)中各分子位置計算系統(tǒng)的勢能,按照經(jīng)典力學(xué)Gear所提出的一種利用數(shù)值解的方法,稱為校正預(yù)測法系統(tǒng)中任一原子i所受的力為勢能的梯度( predictor- corrector method)。時間t+時的位置、速度等可由時間t的泰勒展開式預(yù)測得到(1)r"(+6)=()+8n()+12()+18b(+由牛頓第二定律可得i原子的加速度為F(+8)=()+n()+182b()(t+)=a()+b(1)將牛頓運動定律方程式對時間積分,可預(yù)測i原子經(jīng)過時間t后的速度與位置:b(+)=6()+式中的(、a()、b()為()的1次、2次、3次微分。式(7)dt所產(chǎn)生的速度、加速度不是由牛頓運動方程解得的,所以并非(3)完全正確。可由所預(yù)測的位置產(chǎn)(+)計算所受的力及正確的加速度d(+)。設(shè)正確的加速度與預(yù)測的加速度之間的誤差為r=r +p t+△a(+8)=a(+8)-a(+8)(8)式中,及"分別是粒子i的位置與速度,上標(biāo)“0”為各可得各量的校正式為:物理量的初始值[1]r(+8)=r"(+8)+c△a(+2牛頓運動方程的數(shù)值解法為了得到原子的運動軌跡,必須解式(3)的牛頓運動方(+8)=(+6)+c△(+)(9)程,可采用有限差分法。有限差分法的基本思想就是將積分分a(+8)=a(+6)+c,△(+8)成很多小步,每一小步的時間固定為c。常用的有以下兒種算法:① Verlet算法;② Velocity- Verlet算法:③leapb(+)=b(+)+c3△a(+)frog算法(蛙跳算法);④ Beeman算法:⑤Gear算法式中,C0、C1c2、C3均為常數(shù)。以上為Gear的一次預(yù)測1 eap-frog算法和Gear算法由于使用簡便,準(zhǔn)確性及穩(wěn)定性校正法,也可將此計算推展至更高次的校正高,節(jié)省儲存空間等作者: photon優(yōu)點,已被廣泛采用3勢函數(shù)2.1leap-frog算法勢函數(shù)表明了原子間的相互作用。針對不同的計算物質(zhì)Leap-frog算法速度與位置的數(shù)學(xué)式為:不同的模擬目的,勢函數(shù)有不同的形式。計算結(jié)果的可靠性與勢函數(shù)密切相關(guān)。在分子動力學(xué)發(fā)展初期,主要采用對勢。隨(+-=v(r--S)+sat著模擬體系的復(fù)雜性,逐漸出現(xiàn)了多體勢,以彌補對勢的不r(t+)=r()+Cp(t+-)3.1對勢為了執(zhí)行1eap-frog算法,必須首先由t-0.5時刻的速主要是 Lennard- Jones勢(L-J勢),又叫12-6勢能,它度與t時刻的加速度計算出(t+8),然后由方程的數(shù)學(xué)表達式是(r+)=(+16(+)+5b()-16(t-06)(5)U(r)=4中國煤化工(10)計算出位置r(+6t)。時間為t時的速度可由式(6)算出,即式中,r為CNMHG能參數(shù)能中,r-2項是排斥項,r項是吸引項。當(dāng)r很大時,L-J勢能p()=(+)+p(-O)(6)趨近于零,表示當(dāng)原子對相距很遠(yuǎn)時,彼此之間已經(jīng)沒有非鍵SILICON【高新技術(shù)產(chǎn)業(yè)發(fā)展】 VALLE結(jié)作用了。這個形式的勢函數(shù)表達的作用力較弱,適合描述惰定義為性氣體的固體和液體[2]多體勢Daw和 Baskes首次提出了嵌入原子勢(EAM),其體系的總式中,P為系統(tǒng)密度,dN為距離參考分子中心由r→r+ah勢能可以表示為:球殼內(nèi)的平均粒子數(shù)。徑向分布函數(shù)可以理解為所模擬體系的區(qū)域密度與平均密U=∑F(p)+∑()(11)度的比[1]。當(dāng)r值小時,距離參考分子很近,區(qū)域密度不同于平均密度;當(dāng)r值大時,距離參考分子很遠(yuǎn),區(qū)域密度與系統(tǒng)的式中,第一項F是嵌入能,表示原子核鑲嵌在電子云背平均密度相同,徑向分布函數(shù)接近于1景中的嵌入能;第二項是對勢能,表示位于晶格點陣上的原子2配位數(shù)核之間的相互作用。多體勢大都用來模擬金屬的微觀性質(zhì)配位數(shù)( coordination number)是指某個原子的最近鄰[2]。后來 Baskes等人修正了嵌入原子法,提出MEA勢,可以原子個數(shù),它顯示了該原子周圍的原子分布的密度大小,是微描述共價鍵材料觀結(jié)構(gòu)分析的輔助物理量[4]。其定義為:4周期性邊界條件與最近鏡像(13)執(zhí)行分子動力計算時,通常選取一定數(shù)目N的分子或原n=4mJ”g()2bh子,將其放在一個立方體的盒子中,該盒子即為模擬系統(tǒng)。而式中,P是對應(yīng)點的數(shù)密度,積分的上下限由徑向分布函計算機最多只能模擬兒百到幾千個粒子的系統(tǒng),為了能夠用較數(shù)g(r)得到,當(dāng)萬=0,n=rm時,得到第一配位圈的配位少的粒子數(shù)目來模擬真實的宏觀體系,引入三維周期性邊界條數(shù)件。在三維體系中,每一個單胞會被其它26個單胞所包圍。當(dāng)6.3擴散系數(shù)模擬的單胞中一個粒子由于力的作用離開這個單胞的時候,就分子動力學(xué)模擬中,采用兩種方法來計算體系分子的擴散會有另一個和它對應(yīng)的粒子運動到這個單胞中來,這樣,模擬系數(shù)D,分別是 Einstein法和reen-Kubo法。 Einstein法是利的整個體系的粒子數(shù)就保持不變,密度也不變,符合實際的要對均方位移求斜率來求擴散系數(shù),計算公式為求[3]。D=imaF()-7(0)(14)計算系統(tǒng)中分子間的作用力時,采取最近鏡像的方法。計算分子A與B的作用力,是取與分子A和其最近的距離鏡像分子Green-Kubo法是通過對速度自相關(guān)函數(shù)("ACF)的積分獲B,而非計算系統(tǒng)中的分子A與分子B。因為在計算中利用了最近得擴散系數(shù),計算公式為鏡像的概念,因此必須采用截斷半徑的方法計算非鍵結(jié)的遠(yuǎn)程(15)作用力。定義為截斷半徑,則當(dāng)>時,勢能值趨近于零表示分子間的范德瓦耳斯作用力可忽略不計。截斷半徑最大不其中,D表示粒子i的擴散系數(shù),r()、r(0)分別表示粒能超過盒子邊長的一半,即r