第三章_經(jīng)典分子動(dòng)力學(xué)方法

1、第三章第三章 經(jīng)典分子動(dòng)力學(xué)方法經(jīng)典分子動(dòng)力學(xué)方法3.1 引言引言 分子動(dòng)力學(xué) (Molecular Dynamics, 簡(jiǎn)寫為MD) 方法是確定性模擬方法，這方法是按該體系內(nèi)部的內(nèi)稟動(dòng)力學(xué)規(guī)律來(lái)計(jì)算確定位形的轉(zhuǎn)變。 首先需要建立一組分子的運(yùn)動(dòng)方程，然后通過(guò)直接對(duì)系統(tǒng)中的每一個(gè)原子/分子運(yùn)動(dòng)方程進(jìn)行數(shù)值求解，得到每個(gè)時(shí)刻每個(gè)原子/分子的坐標(biāo)與動(dòng)量(速度)，即在相空間的運(yùn)動(dòng)軌跡，再利用統(tǒng)計(jì)方法得到多體系統(tǒng)的靜態(tài)和動(dòng)態(tài)特性，從而得到系統(tǒng)的宏觀性質(zhì)。 在MD方法的處理過(guò)程中，方程組的建立是通過(guò)對(duì)物理體系的微觀數(shù)學(xué)描述給出的，在這個(gè)微觀的物理體系中，每個(gè)原子/分子都各自服從經(jīng)典的牛頓力學(xué)定律。 MD方

2、法是實(shí)現(xiàn)玻爾茲曼的統(tǒng)計(jì)力學(xué)途徑，可以處理與時(shí)間有關(guān)的過(guò)程，因而可以處理非平衡態(tài)問題，但是該方法的計(jì)算機(jī)程序較復(fù)雜，計(jì)算量大，占內(nèi)存也多。MD方法的發(fā)展史方法的發(fā)展史 MD方法是20世紀(jì)50年代后期由B.J Alder和T.E. Wainwright創(chuàng)造發(fā)展的。他們?cè)?957年利用MD方法，發(fā)現(xiàn)了早在1939年根據(jù)統(tǒng)計(jì)力學(xué)預(yù)言的“剛性球組成的集合系統(tǒng)會(huì)發(fā)生由其液相到結(jié)晶相的相轉(zhuǎn)變”。 20世紀(jì)70年代，產(chǎn)生了剛性體系的動(dòng)力學(xué)方法被應(yīng)用于水和氮等分子性溶液體系的處理，取得了成功。1972年，A.W. Less和S.F. Edwards等人發(fā)展了該方法，并擴(kuò)展到了存在速度梯度(即處于非平衡狀態(tài))的系

3、統(tǒng)。 之后，此方法被M.J. Gillan等人推廣到了具有溫度梯度的非平衡系統(tǒng)，從而構(gòu)造并形成了非平衡MD方法體系。MD方法的發(fā)展史方法的發(fā)展史 MD方法真正作為材料科學(xué)領(lǐng)域的一個(gè)重要研究方法，開始于恒壓MD方法(1980)和恒溫MD方法(1984)的建立及在應(yīng)用方面的成功。 1985年人們又提出了將電子論和分子動(dòng)力學(xué)方法有機(jī)統(tǒng)一起來(lái)的所謂CarParrinello方法，即第一性原理MD方法。它不僅可以處理半導(dǎo)體和金屬的問題，同時(shí)還可應(yīng)用于處理有機(jī)物和化學(xué)反應(yīng)。 1991年有人進(jìn)一步提出了巨正則系綜MD方法，從而又可適用于吸附問題的處理等，該方法還在進(jìn)一步發(fā)展之中。 分子動(dòng)力學(xué)方法的主要發(fā)展可

4、見表3.1。年代創(chuàng)立者創(chuàng)造內(nèi)容工作（MD分類名稱）1957B. J Alder & T. E. Wainwright剛性球MD方法1963A. Rahman質(zhì)點(diǎn)系MD方法1971Rahman &F. H. Stillinger剛性系統(tǒng)MD方法1972W. Lees &S. F. Edwards平衡系統(tǒng)MD方法（存在速度梯度）1977J. P. Rychaert et al.約束系統(tǒng)MD方法1980Andersen, Parrinello & Rahman恒壓恒壓MD方法方法1983N. J. Gillan &M. Dixon非平衡MD方法（存在溫度梯度）

5、1984 S. Nos恒溫恒溫MD方法方法1985R. Car &M. Parrinello第一性原理第一性原理MD方法方法(Car-Parrinello方法方法)1991Cagin & Pettitt巨正則系統(tǒng)MD方法表3.1 MD方法的里程碑工作方法的里程碑工作 3.2 MD方法計(jì)算初步方法計(jì)算初步 在計(jì)算機(jī)出現(xiàn)以前，作為根據(jù)原子間相互作用力等微觀信息了解多原子或分子團(tuán)的結(jié)構(gòu)、性質(zhì)的方法，所采用的是基于統(tǒng)計(jì)理論的數(shù)學(xué)解析法。然而，原子間相互作用力稍微復(fù)雜一些，不用說(shuō)求解統(tǒng)計(jì)理論嚴(yán)格方程解，就是進(jìn)行數(shù)值求解也是一件很困難的事。 MD方法就是數(shù)值求解多體系統(tǒng)的確定性運(yùn)動(dòng)方程，并

6、根據(jù)對(duì)所求結(jié)果進(jìn)行統(tǒng)計(jì)處理，決定粒子的軌跡，從而給出物性預(yù)測(cè)和微觀結(jié)構(gòu)信息的一種模擬方法。內(nèi)能比熱容運(yùn)動(dòng)方程溫度、壓力相互作用原子位置坐標(biāo)3維結(jié)構(gòu)原子坐標(biāo)、速度原子運(yùn)動(dòng)熱力學(xué)性質(zhì)動(dòng)力學(xué)性質(zhì)光學(xué)性質(zhì)（輸出信息）（二次信息）擴(kuò)散系數(shù)粘滯系數(shù)電導(dǎo)率紅外吸收?qǐng)D圖3-1 MD方法信息輸入輸出信息方框圖方法信息輸入輸出信息方框圖 MD這種方法并不嚴(yán)格。因此，必須根據(jù)情況，檢驗(yàn)改變所模擬的基本單元尺寸所得結(jié)果是否會(huì)改變，直到所得結(jié)果不隨基本單元尺寸變化而變化。通常這樣的處理在很多情況下是有效的。 對(duì)于基本單元內(nèi)的原子、分子運(yùn)動(dòng)方程，使用什么樣的形式合適，要具體問題具體分析。若是考慮具有確定的粒子數(shù)N，體積V

7、和能量E的NEV系綜(稱為微正則系綜，Micro-Canonical Ensemble)，則其運(yùn)動(dòng)方程可以表達(dá)成式(3-2-1)所示的普通牛頓方程的形式22iiid rmFdt (3-2-1)式中mi為所考察的原子質(zhì)量，ri為原子的位置坐標(biāo)，F(xiàn)i為作用在原子上的原子相互作用的合力，它由下式給出1NiiijjF (3-2-2)其中，ij是原子和原子j之間的勢(shì)函數(shù)(有時(shí)亦稱為力場(chǎng)) 例如，由氬原子等組成的稀薄氣體，其勢(shì)函數(shù)可采用LennardJones勢(shì)，1264ijijijijrr (3-2-3)式中，r 是原子間距， 是結(jié)合強(qiáng)度參數(shù)， 是表示原子半徑的參數(shù)。ijij在t時(shí)間內(nèi)，對(duì)系統(tǒng)內(nèi)的所有粒

8、子解運(yùn)動(dòng)方程111()()()()ininininnnr tr tvtvtttt 1nmaxttYesNo啟動(dòng)計(jì)算設(shè)定坐標(biāo)、速度初始值計(jì)算作用在原子上的力Fi計(jì)算要求的物理量，將數(shù)據(jù)寫入軌跡文件 ttmax輸出計(jì)算結(jié)果，并結(jié)束計(jì)算對(duì)(3-2-1)可用數(shù)值積分法求解，其數(shù)據(jù)處理流程圖見圖3-2 圖圖3-2 MD數(shù)據(jù)處理流程圖數(shù)據(jù)處理流程圖MD方法NEV能量恒定NTV恒溫NHP恒壓NTP恒溫恒壓恒溫恒壓VT巨正則系VL恒化學(xué)勢(shì)彈性力學(xué)（原子分子）質(zhì)點(diǎn)力學(xué)質(zhì)點(diǎn)力學(xué)（原子分子）剛體力學(xué)（分子）約束力（分子和晶體）系）動(dòng)力學(xué)模型目標(biāo)系統(tǒng)團(tuán)簇塊體材料塊體材料表面界面MBE/CVDBerrele法Green

9、法多重時(shí)間刻度法數(shù)值積分法邊界條件到目前為止已經(jīng)確立的MD方法的方法的主要技術(shù)體系 統(tǒng)計(jì)系綜統(tǒng)計(jì)系綜 系綜是一個(gè)巨大的系統(tǒng)，由組成、性質(zhì)、尺寸和形狀完全一樣的全同體系構(gòu)成數(shù)目極多的系統(tǒng)的集合。不同的系綜，MD方法的基本方程有所不同。 目前除微正則系綜（NEV系綜）外，已完成了正則系綜（NTV系綜），等溫等壓系綜（NTP系綜），等壓等焓系綜（NHP系綜），巨正則系綜（VL系綜），恒定化學(xué)勢(shì)系綜（VT系綜）等五個(gè)系綜的MD方法的基本方程的確立。 現(xiàn)在已經(jīng)能夠處理許多體系，例如： 孤立宏觀團(tuán)簇的模擬(用NEV或NTV系綜) 固體的結(jié)構(gòu)相變，玻璃轉(zhuǎn)變，晶化過(guò)程的模擬(用NTP系綜) 固體(晶體)表面的

10、原子、分子吸附現(xiàn)象的模擬 (用VT系綜)力學(xué)條件力學(xué)條件 已建立了彈性力學(xué)、質(zhì)點(diǎn)力學(xué)、剛體力學(xué)、約束力學(xué)等不同力學(xué)條件下的四種體系的MD方法。 彈性力學(xué)方法是將所考察的原子分子看作剛性球來(lái)處理，建立完全彈性碰撞方程，借以求解出原子、分子的運(yùn)動(dòng)規(guī)律。這種處理可以在液晶的模擬中使用。 質(zhì)點(diǎn)力學(xué)模型質(zhì)點(diǎn)力學(xué)模型是將原子、分子作為質(zhì)點(diǎn)處理，粒子間的相互作用力采用坐標(biāo)的連續(xù)函數(shù)。這種力學(xué)體系的應(yīng)用對(duì)象非常多，可以用于處理陶瓷、金屬、半導(dǎo)體等無(wú)機(jī)化合物材料以及有機(jī)高分子、生物大分子等幾乎所有的材料。力學(xué)條件力學(xué)條件 剛體力學(xué)方法是把分子作為剛體處理，建立對(duì)于剛體的歐拉方程和對(duì)于質(zhì)點(diǎn)系的牛頓方程，聯(lián)立求解所

11、研究對(duì)象的運(yùn)動(dòng)問題，以前主要是在處理像水和四氯化碳那樣的低分子量體系，現(xiàn)已用于研究晶體的相變。 約束力學(xué)是凍結(jié)粒子體系的一部分自由度，進(jìn)而求解因此而生成約束條件下的質(zhì)點(diǎn)力學(xué)運(yùn)動(dòng)方程。對(duì)于有機(jī)分子來(lái)講，因?yàn)殒I長(zhǎng)、鍵角的振動(dòng)變化對(duì)體系的結(jié)構(gòu)影響不大，將這些自由度凍結(jié)是合適的。另外也有采用固定晶體結(jié)構(gòu)的考慮。邊界條件問題邊界條件問題 在處理原子、分子的聚合體問題時(shí)，就MD方法而言，能處理的原子(分子)數(shù)目要受到計(jì)算機(jī)運(yùn)行速度和能力的限制。目前報(bào)導(dǎo)的最好水平是能處理109量級(jí)的原子數(shù)目。這與現(xiàn)實(shí)物質(zhì)含有1023個(gè)原子或分子的差距還很大，導(dǎo)致模擬系統(tǒng)原子數(shù)少于真實(shí)系統(tǒng)的所謂“尺寸效應(yīng)”的問題。 為了減小

12、“尺寸效應(yīng)”而又不至于使計(jì)算工作量過(guò)大，對(duì)于平衡態(tài)MD模擬采用 “周期性邊界條件”。周期性邊界條件 圖3-4假設(shè)現(xiàn)實(shí)物質(zhì)中的一部分原子(通常為102105個(gè))，被取出配置在所謂基本單元的箱中，由于基本單元周圍的原子、分子變成表面，從而不同于本來(lái)要處理的體狀態(tài)。為了防止這種情況，在基本單元周圍配置其復(fù)制品(圖3-4)。周期性邊界條件 對(duì)于在基本單元周圍邊界的原子、分子所受到的作用力，不僅有來(lái)自基本單元之內(nèi)部的原子或分子的作用力的貢獻(xiàn)，還要考慮來(lái)自其近鄰復(fù)制單元的原子或分子的作用力的貢獻(xiàn)。這樣，把來(lái)自距粒子某一距離（截?cái)嗑嚯x）內(nèi)的j粒子的貢獻(xiàn)（這里與粒子j是否處在基本單元內(nèi)或復(fù)制單元內(nèi)無(wú)關(guān)）求和。

13、給出力的方法稱為最小鏡像距離法（Minimum Image Distance）。周期性邊界條件如果在基本單元內(nèi)的原子的位置為, 周期邊界條件會(huì)產(chǎn)生該原子的鏡像。它的位置在 這兒a、b、c是基本單元的三個(gè)邊長(zhǎng)，l，m，n是整數(shù)，取值范圍為- 到+。 在基本單元內(nèi)的原子不僅與在本單元內(nèi)的其他原子有相互作用，還與在相鄰單元內(nèi)的鏡像原子有作用。imageiirrlambnc (3-2-4)基本單元大小的選擇基本單元大小的選擇 基本單元的大小必須大于2Rcut（Rcut是相互作用勢(shì)的截?cái)嗑嚯x）或Rcut1012 m2/s時(shí)適用。 2211NiiiMSDr tr tr oN (3-3-10)(3-3-9)

14、GreenKubo formula 擴(kuò)散系數(shù)D也可用速度自相關(guān)函數(shù)通過(guò)下面的GreenKubo 公式表示 （3311） GreenKubo公式和Einstein 關(guān)系是等價(jià)的。 101133NiiiDv tv odtvtvo dtNPartial diffusion 在多種原子組成的體系中（例如AB2），可以對(duì)某一種原子求和，得到Partial diffusion，從而得到原子A或原子B的擴(kuò)散系數(shù)DA、DB。非均勻體系中的擴(kuò)散非均勻體系中的擴(kuò)散 對(duì)一非均勻體系，例如一塊材料，有體和晶界部分?？煞謩e對(duì)某一類原子做平均，如對(duì)晶界原子，或?qū)w中的原子，得到不同部分原子不同的擴(kuò)散。與溫度有關(guān)的擴(kuò)散與溫

15、度有關(guān)的擴(kuò)散 在一個(gè)系統(tǒng)里，若原子的擴(kuò)散是通過(guò)跳躍進(jìn)行的，可將擴(kuò)散用Arrhenius關(guān)系 來(lái)表達(dá)，式中是prefactor，是原子遷移率或激發(fā)能(activation energy)。對(duì)（3312）兩邊求對(duì)數(shù)有0expaBEDDk T(3-3-12)TkEDDBa1lnln0(3-3-13)與溫度有關(guān)的擴(kuò)散與溫度有關(guān)的擴(kuò)散 作lnD1/T圖(圖3.8)，從圖中直線的斜率可以得到Ea，從截距可以得到D0。同理也可以得到表面原子、晶界原子的擴(kuò)散系數(shù)(圖3.9)。圖3.8圖3.9TkEDDBa1lnln0MD模擬材料的溶解模擬材料的溶解 固態(tài)到液態(tài)的相變（溶解）的過(guò)程發(fā)生在幾ps時(shí)間，在MD模擬中可通過(guò)以下手段知道是否發(fā)生了