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