第三章

經(jīng)典分子動(dòng)力學(xué)方法1精選ppt課件20213.1引言分子動(dòng)力學(xué)(MolecularDynamics,簡(jiǎn)寫為MD)方法是確定性模擬方法，這方法是按該體系內(nèi)部的內(nèi)稟動(dòng)力學(xué)規(guī)律來計(jì)算確定位形的轉(zhuǎn)變。首先需要建立一組分子的運(yùn)動(dòng)方程，然后通過直接對(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方法的處理過程中，方程組的建立是通過對(duì)物理體系的微觀數(shù)學(xué)描述給出的，在這個(gè)微觀的物理體系中，每個(gè)原子/分子都各自服從經(jīng)典的牛頓力學(xué)定律。MD方法是實(shí)現(xiàn)玻爾茲曼的統(tǒng)計(jì)力學(xué)途徑，可以處理與時(shí)間有關(guān)的過程，因而可以處理非平衡態(tài)問題，但是該方法的計(jì)算機(jī)程序較復(fù)雜，計(jì)算量大，占內(nèi)存也多。2精選ppt課件2021MD方法的發(fā)展史MD方法是20世紀(jì)50年代后期由B.JAlder和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))的系統(tǒng)。之后，此方法被M.J.Gillan等人推廣到了具有溫度梯度的非平衡系統(tǒng)，從而構(gòu)造并形成了非平衡MD方法體系。3精選ppt課件2021MD方法的發(fā)展史MD方法真正作為材料科學(xué)領(lǐng)域的一個(gè)重要研究方法，開始于恒壓MD方法(1980)和恒溫MD方法(1984)的建立及在應(yīng)用方面的成功。1985年人們又提出了將電子論和分子動(dòng)力學(xué)方法有機(jī)統(tǒng)一起來的所謂Car－Parrinello方法，即第一性原理MD方法。它不僅可以處理半導(dǎo)體和金屬的問題，同時(shí)還可應(yīng)用于處理有機(jī)物和化學(xué)反應(yīng)。1991年有人進(jìn)一步提出了巨正則系綜MD方法，從而又可適用于吸附問題的處理等，該方法還在進(jìn)一步發(fā)展之中。分子動(dòng)力學(xué)方法的主要發(fā)展可見表3.1。4精選ppt課件2021年代創(chuàng)立者創(chuàng)造內(nèi)容工作（MD分類名稱）1957B.JAlder&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.Rychaertetal.約束系統(tǒng)MD方法1980Andersen,Parrinello&Rahman恒壓MD方法1983N.J.Gillan&M.Dixon非平衡MD方法（存在溫度梯度）1984S.Nosé恒溫MD方法1985R.Car&M.Parrinello第一性原理MD方法(Car-Parrinello方法)1991Cagin&Pettitt巨正則系統(tǒng)MD方法表3.1MD方法的里程碑工作

5精選ppt課件20213.2MD方法計(jì)算初步在計(jì)算機(jī)出現(xiàn)以前，作為根據(jù)原子間相互作用力等微觀信息了解多原子或分子團(tuán)的結(jié)構(gòu)、性質(zhì)的方法，所采用的是基于統(tǒng)計(jì)理論的數(shù)學(xué)解析法。然而，原子間相互作用力稍微復(fù)雜一些，不用說求解統(tǒng)計(jì)理論嚴(yán)格方程解，就是進(jìn)行數(shù)值求解也是一件很困難的事。MD方法就是數(shù)值求解多體系統(tǒng)的確定性運(yùn)動(dòng)方程，并根據(jù)對(duì)所求結(jié)果進(jìn)行統(tǒng)計(jì)處理，決定粒子的軌跡，從而給出物性預(yù)測(cè)和微觀結(jié)構(gòu)信息的一種模擬方法。6精選ppt課件2021內(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)D3-1MD方法信息輸入輸出信息方框圖7精選ppt課件2021MD這種方法并不嚴(yán)格。因此，必須根據(jù)情況，檢驗(yàn)改變所模擬的基本單元尺寸所得結(jié)果是否會(huì)改變，直到所得結(jié)果不隨基本單元尺寸變化而變化。通常這樣的處理在很多情況下是有效的。8精選ppt課件2021

對(duì)于基本單元內(nèi)的原子、分子運(yùn)動(dòng)方程，使用什么樣的形式合適，要具體問題具體分析。若是考慮具有確定的粒子數(shù)N，體積V和能量E的NEV系綜(稱為微正則系綜，Micro-CanonicalEnsemble)，則其運(yùn)動(dòng)方程可以表達(dá)成式(3-2-1)所示的普通牛頓方程的形式

(3-2-1)式中mi為所考察的原子質(zhì)量，ri為原子的位置坐標(biāo)，F(xiàn)i為作用在原子上的原子相互作用的合力，它由下式給出

(3-2-2)其中，Φij是原子和原子j之間的勢(shì)函數(shù)(有時(shí)亦稱為力場(chǎng))9精選ppt課件2021例如，由氬原子等組成的稀薄氣體，其勢(shì)函數(shù)可采用Lennard－Jones勢(shì)，

(3-2-3)式中，r是原子間距，是結(jié)合強(qiáng)度參數(shù)，是表示原子半徑的參數(shù)。10精選ppt課件2021在Δt時(shí)間內(nèi)，對(duì)系統(tǒng)內(nèi)的所有粒子解運(yùn)動(dòng)方程YesNo啟動(dòng)計(jì)算設(shè)定坐標(biāo)、速度初始值計(jì)算作用在原子上的力Fi計(jì)算要求的物理量，將數(shù)據(jù)寫入軌跡文件t>tmax輸出計(jì)算結(jié)果，并結(jié)束計(jì)算對(duì)(3-2-1)可用數(shù)值積分法求解，其數(shù)據(jù)處理流程圖見圖3-2圖3-2MD數(shù)據(jù)處理流程圖11精選ppt課件2021MD方法NEV能量恒定NTV恒溫NHP恒壓NTP恒溫恒壓

VT巨正則系

VL恒化學(xué)勢(shì)彈性力學(xué)（原子分子）質(zhì)點(diǎn)力學(xué)（原子分子）剛體力學(xué)（分子）約束力（分子和晶體）系）動(dòng)力學(xué)模型目標(biāo)系統(tǒng)團(tuán)簇塊體材料表面界面MBE/CVDBerrele法Green法多重時(shí)間刻度法數(shù)值積分法邊界條件到目前為止已經(jīng)確立的MD方法的主要技術(shù)體系12精選ppt課件2021①統(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)變，晶化過程的模擬(用NTP系綜)固體(晶體)表面的原子、分子吸附現(xiàn)象的模擬(用μVT系綜)13精選ppt課件2021②力學(xué)條件

已建立了彈性力學(xué)、質(zhì)點(diǎn)力學(xué)、剛體力學(xué)、約束力學(xué)等不同力學(xué)條件下的四種體系的MD方法。彈性力學(xué)方法是將所考察的原子分子看作剛性球來處理，建立完全彈性碰撞方程，借以求解出原子、分子的運(yùn)動(dòng)規(guī)律。這種處理可以在液晶的模擬中使用。質(zhì)點(diǎn)力學(xué)模型是將原子、分子作為質(zhì)點(diǎn)處理，粒子間的相互作用力采用坐標(biāo)的連續(xù)函數(shù)。這種力學(xué)體系的應(yīng)用對(duì)象非常多，可以用于處理陶瓷、金屬、半導(dǎo)體等無機(jī)化合物材料以及有機(jī)高分子、生物大分子等幾乎所有的材料。14精選ppt課件2021②力學(xué)條件剛體力學(xué)方法是把分子作為剛體處理，建立對(duì)于剛體的歐拉方程和對(duì)于質(zhì)點(diǎn)系的牛頓方程，聯(lián)立求解所研究對(duì)象的運(yùn)動(dòng)問題，以前主要是在處理像水和四氯化碳那樣的低分子量體系，現(xiàn)已用于研究晶體的相變。約束力學(xué)是凍結(jié)粒子體系的一部分自由度，進(jìn)而求解因此而生成約束條件下的質(zhì)點(diǎn)力學(xué)運(yùn)動(dòng)方程。對(duì)于有機(jī)分子來講，因?yàn)殒I長(zhǎng)、鍵角的振動(dòng)變化對(duì)體系的結(jié)構(gòu)影響不大，將這些自由度凍結(jié)是合適的。另外也有采用固定晶體結(jié)構(gòu)的考慮。15精選ppt課件2021③邊界條件問題在處理原子、分子的聚合體問題時(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)”的問題。為了減小“尺寸效應(yīng)”而又不至于使計(jì)算工作量過大，對(duì)于平衡態(tài)MD模擬采用“周期性邊界條件”。16精選ppt課件2021周期性邊界條件圖3-4假設(shè)現(xiàn)實(shí)物質(zhì)中的一部分原子(通常為102－105個(gè))，被取出配置在所謂基本單元的箱中，由于基本單元周圍的原子、分子變成表面，從而不同于本來要處理的體狀態(tài)。為了防止這種情況，在基本單元周圍配置其復(fù)制品(圖3-4)。17精選ppt課件2021周期性邊界條件對(duì)于在基本單元周圍邊界的原子、分子所受到的作用力，不僅有來自基本單元之內(nèi)部的原子或分子的作用力的貢獻(xiàn)，還要考慮來自其近鄰復(fù)制單元的原子或分子的作用力的貢獻(xiàn)。這樣，把來自距粒子某一距離（截?cái)嗑嚯x）內(nèi)的j粒子的貢獻(xiàn)（這里與粒子j是否處在基本單元內(nèi)或復(fù)制單元內(nèi)無關(guān)）求和。給出力的方法稱為最小鏡像距離法（MinimumImageDistance）。18精選ppt課件2021周期性邊界條件如果在基本單元內(nèi)的原子的位置為,周期邊界條件會(huì)產(chǎn)生該原子的鏡像。它的位置在

這兒a、b、c是基本單元的三個(gè)邊長(zhǎng)，l，m，n是整數(shù)，取值范圍為-∞到+∞。

在基本單元內(nèi)的原子不僅與在本單元內(nèi)的其他原子有相互作用，還與在相鄰單元內(nèi)的鏡像原子有作用。(3-2-4)19精選ppt課件2021基本單元大小的選擇

基本單元的大小必須大于2Rcut（Rcut是相互作用勢(shì)的截?cái)嗑嚯x）或Rcut<1/2基本單元的大小。這保證了任何原子只與原子的一個(gè)鏡像有相互作用，不與自己的鏡像作用。這個(gè)條件稱為“minimumimagecriterion”在我們所研究的體系內(nèi)的任何結(jié)構(gòu)特性的特征尺寸或任何重要的效應(yīng)的特征長(zhǎng)度必須小于基本單元的大小。為了檢驗(yàn)不同基本單元大小是否會(huì)引入“人為效應(yīng)”，必須用不同的基本單元尺寸做計(jì)算，若結(jié)果能收斂，則尺寸選擇是合適的。20精選ppt課件2021邊界條件邊界條件的問題比較復(fù)雜，因所考慮的具體情況而定，概括地說，作為處理原子、分子團(tuán)簇的邊界條件可進(jìn)行推廣。例如，適用于塊體狀態(tài)的三維周期邊界條件，可擴(kuò)展處理表面重構(gòu)和表面吸附、單層膜相變的二維周期邊界條件，以及處理異質(zhì)晶體界面的邊界系統(tǒng)等，同時(shí)還開發(fā)了模擬非平衡的方法。諸如針對(duì)分子束外延（MBE）的邊界條件，人們進(jìn)行了很多研究。以二維邊界條件為例。只在x-y平面配置基本單元的復(fù)制品，使用周期性邊界條件，同時(shí)在z方向不賦予周期性邊界條件，固定兩端的數(shù)層原子。由于采用了這樣的人工邊界條件，所以在z軸方向上的原子層數(shù)要有適當(dāng)?shù)臄?shù)目，一般考慮的標(biāo)準(zhǔn)線度為4－5納米。采用二維邊界條件，就可以用MD在原子尺度上研究界面的結(jié)構(gòu)。21精選ppt課件2021④數(shù)值積分法MD方法的基本方程是線性或非線性的二階常微分方程。對(duì)此人們研究了許多求解方法。例如，貝魯勒（Berreele）法，阿達(dá)姆斯（Adams）法，龍格-庫(kù)塔（Runge-Kuta）法，追趕法等。最近，M.Tuckerman和B.J.Berne提出了所謂多重時(shí)間寬度計(jì)算方法。即在質(zhì)量相差很大的多原子體系，力學(xué)常數(shù)有很大不同的體系以及含有長(zhǎng)、短程力的體系等情況下，必須結(jié)合具有最小振動(dòng)周期的自由度選擇積分的時(shí)間標(biāo)度（）.這方法縮短了計(jì)算用的時(shí)間，有望成為MD方法中有前途的數(shù)值積分法。通常～10-15s=1fs22精選ppt課件20213.3MD模擬的數(shù)據(jù)分析

MD模擬給出的直接數(shù)據(jù)是坐標(biāo)和速度，對(duì)它們進(jìn)行加工，可得到各種所要的信息。這些信息大致可分為3類：·所模擬系統(tǒng)的平衡態(tài)性質(zhì)·系統(tǒng)處于亞穩(wěn)態(tài)時(shí)的結(jié)構(gòu)與性質(zhì)·系統(tǒng)在遠(yuǎn)離平衡態(tài)時(shí)的動(dòng)力學(xué)過程23精選ppt課件2021平衡態(tài)一個(gè)物理量A隨時(shí)間變化而趨于平衡態(tài)的值A(chǔ)0的過程可用下式表示

其中

是馳豫時(shí)間。對(duì)于小的

，我們只需等待系統(tǒng)達(dá)到平衡態(tài)后，再開始收集系統(tǒng)在平衡態(tài)時(shí)的參量。對(duì)很大的

，要達(dá)平衡態(tài)，需要模擬很長(zhǎng)的時(shí)間，這時(shí)MD方法是不合適的。對(duì)中等大小的

，如果我們不能直接得到，但可以估計(jì)A0。在許多情況下，我們并不需要達(dá)到平衡。如果我們的目的是研究非平衡的過程。(3-3-1)24精選ppt課件2021在MD模擬中可觀察量的測(cè)量可把可觀察量表示為一個(gè)位置和速度的函數(shù)

對(duì)上式做時(shí)間平均，得

我們更關(guān)心的是整個(gè)系統(tǒng)的平均效應(yīng)。(3-3-2)(3-3-3)25精選ppt課件2021為什么要計(jì)算能量？驗(yàn)證體系是否總能量守恒若有能量從動(dòng)能轉(zhuǎn)移到勢(shì)能，表示體系發(fā)生了相變?nèi)粼贓～T曲線上有跳躍，表示有一個(gè)一階相變。26精選ppt課件2021能量計(jì)算每個(gè)原子的勢(shì)能為（3－3－4）每個(gè)原子的動(dòng)能為（3－3－5）每個(gè)原子的總能量（3－3－6）27精選ppt課件2021一般來說，對(duì)于給定能量的精確初始條件是無法知道的，為了把系統(tǒng)調(diào)節(jié)到給定的能量，先給出一個(gè)合理的初始條件，然后對(duì)能量進(jìn)行增減，直至系統(tǒng)達(dá)到所要達(dá)到的狀態(tài)為止。能量的調(diào)整一般是通過對(duì)速度進(jìn)行特別的標(biāo)度(scaling)來實(shí)現(xiàn)的。這種標(biāo)度可以使系統(tǒng)的速度發(fā)生很大的變化。為了消除可能帶來的效應(yīng)，必須給系統(tǒng)足夠的時(shí)間以再次建立平衡。其具體步驟為：1.解運(yùn)動(dòng)方程，給出一定時(shí)間步的結(jié)果；2.計(jì)算體系的動(dòng)能和勢(shì)能；3.觀察體系的總能量是否守恒，28精選ppt課件2021如總能量不守恒，則通過調(diào)節(jié)速度來實(shí)現(xiàn)，即將速度乘以一個(gè)標(biāo)定因子

：,總能量為:

g為總自由度數(shù)。標(biāo)度因子近似可表達(dá)成

在MD計(jì)算時(shí)，這個(gè)給系統(tǒng)足夠的時(shí)間以再次建立平衡的時(shí)間稱為scallingtime（一般為幾千到幾十萬個(gè)step），要求在這個(gè)期間內(nèi)，能量守恒。過了這個(gè)時(shí)間，我們才開始采取數(shù)據(jù)做統(tǒng)計(jì)計(jì)算。29精選ppt課件2021溫度T的計(jì)算在平衡時(shí)，速度分布滿足麥克斯韋-波爾茲曼分布，有

這兒是玻爾茲曼常數(shù)，從而得到

可見溫度是原子、分子運(yùn)動(dòng)的統(tǒng)計(jì)平均的結(jié)果(3-3-7)(3-3-8)30精選ppt課件2021原子的遷移率與擴(kuò)散在MD方法中我們用原子移動(dòng)距離的平方的平均值MSD（meansquaredisplacement）來表示原子的遷移率。

我們用Einstein關(guān)系

MSD=A+6Dt+fluctuation（3D系統(tǒng)）Einstein關(guān)系將MSD與擴(kuò)散系數(shù)D相聯(lián)系，從MSD－t曲線的斜率6D可決定擴(kuò)散系數(shù)D。對(duì)3D系統(tǒng)Einstein關(guān)系里Dt前的系數(shù)為6，2D系統(tǒng)該系數(shù)為4,1D系統(tǒng)該系數(shù)為2。這個(gè)關(guān)系只在足夠高的溫度下適用，或當(dāng)D>1012m2/s時(shí)適用。(3-3-10)(3-3-9)31精選ppt課件2021Green－Kuboformula擴(kuò)散系數(shù)D也可用速度自相關(guān)函數(shù)通過下面的Green－Kubo公式表示

（3－3－11）Green－Kubo公式和Einstein關(guān)系是等價(jià)的。32