第三章

經典分子動力學方法1精選ppt課件20213.1引言分子動力學(MolecularDynamics,簡寫為MD)方法是確定性模擬方法，這方法是按該體系內部的內稟動力學規律來計算確定位形的轉變。首先需要建立一組分子的運動方程，然后通過直接對系統中的每一個原子/分子運動方程進行數值求解，得到每個時刻每個原子/分子的坐標與動量(速度)，即在相空間的運動軌跡，再利用統計方法得到多體系統的靜態和動態特性，從而得到系統的宏觀性質。在MD方法的處理過程中，方程組的建立是通過對物理體系的微觀數學描述給出的，在這個微觀的物理體系中，每個原子/分子都各自服從經典的牛頓力學定律。MD方法是實現玻爾茲曼的統計力學途徑，可以處理與時間有關的過程，因而可以處理非平衡態問題，但是該方法的計算機程序較復雜，計算量大，占內存也多。2精選ppt課件2021MD方法的發展史MD方法是20世紀50年代后期由B.JAlder和T.E.Wainwright創造發展的。他們在1957年利用MD方法，發現了早在1939年根據統計力學預言的“剛性球組成的集合系統會發生由其液相到結晶相的相轉變”。20世紀70年代，產生了剛性體系的動力學方法被應用于水和氮等分子性溶液體系的處理，取得了成功。1972年，A.W.Less和S.F.Edwards等人發展了該方法，并擴展到了存在速度梯度(即處于非平衡狀態)的系統。之后，此方法被M.J.Gillan等人推廣到了具有溫度梯度的非平衡系統，從而構造并形成了非平衡MD方法體系。3精選ppt課件2021MD方法的發展史MD方法真正作為材料科學領域的一個重要研究方法，開始于恒壓MD方法(1980)和恒溫MD方法(1984)的建立及在應用方面的成功。1985年人們又提出了將電子論和分子動力學方法有機統一起來的所謂Car－Parrinello方法，即第一性原理MD方法。它不僅可以處理半導體和金屬的問題，同時還可應用于處理有機物和化學反應。1991年有人進一步提出了巨正則系綜MD方法，從而又可適用于吸附問題的處理等，該方法還在進一步發展之中。分子動力學方法的主要發展可見表3.1。4精選ppt課件2021年代創立者創造內容工作（MD分類名稱）1957B.JAlder&T.E.Wainwright剛性球MD方法1963A.Rahman質點系MD方法1971Rahman&F.H.Stillinger剛性系統MD方法1972W.Lees&S.F.Edwards平衡系統MD方法（存在速度梯度）1977J.P.Rychaertetal.約束系統MD方法1980Andersen,Parrinello&Rahman恒壓MD方法1983N.J.Gillan&M.Dixon非平衡MD方法（存在溫度梯度）1984S.Nosé恒溫MD方法1985R.Car&M.Parrinello第一性原理MD方法(Car-Parrinello方法)1991Cagin&Pettitt巨正則系統MD方法表3.1MD方法的里程碑工作

5精選ppt課件20213.2MD方法計算初步在計算機出現以前，作為根據原子間相互作用力等微觀信息了解多原子或分子團的結構、性質的方法，所采用的是基于統計理論的數學解析法。然而，原子間相互作用力稍微復雜一些，不用說求解統計理論嚴格方程解，就是進行數值求解也是一件很困難的事。MD方法就是數值求解多體系統的確定性運動方程，并根據對所求結果進行統計處理，決定粒子的軌跡，從而給出物性預測和微觀結構信息的一種模擬方法。6精選ppt課件2021內能比熱容運動方程溫度、壓力相互作用原子位置坐標3維結構原子坐標、速度原子運動熱力學性質動力學性質光學性質（輸出信息）（二次信息）擴散系數粘滯系數電導率紅外吸收圖3-1MD方法信息輸入輸出信息方框圖7精選ppt課件2021MD這種方法并不嚴格。因此，必須根據情況，檢驗改變所模擬的基本單元尺寸所得結果是否會改變，直到所得結果不隨基本單元尺寸變化而變化。通常這樣的處理在很多情況下是有效的。8精選ppt課件2021

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

(3-2-1)式中mi為所考察的原子質量，ri為原子的位置坐標，Fi為作用在原子上的原子相互作用的合力，它由下式給出

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

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

VT巨正則系

VL恒化學勢彈性力學（原子分子）質點力學（原子分子）剛體力學（分子）約束力（分子和晶體）系）動力學模型目標系統團簇塊體材料表面界面MBE/CVDBerrele法Green法多重時間刻度法數值積分法邊界條件到目前為止已經確立的MD方法的主要技術體系12精選ppt課件2021①統計系綜系綜是一個巨大的系統，由組成、性質、尺寸和形狀完全一樣的全同體系構成數目極多的系統的集合。不同的系綜，MD方法的基本方程有所不同。目前除微正則系綜（NEV系綜）外，已完成了正則系綜（NTV系綜），等溫等壓系綜（NTP系綜），等壓等焓系綜（NHP系綜），巨正則系綜（μVL系綜），恒定化學勢系綜（μVT系綜）等五個系綜的MD方法的基本方程的確立。現在已經能夠處理許多體系，例如：孤立宏觀團簇的模擬(用NEV或NTV系綜)固體的結構相變，玻璃轉變，晶化過程的模擬(用NTP系綜)固體(晶體)表面的原子、分子吸附現象的模擬(用μVT系綜)13精選ppt課件2021②力學條件

已建立了彈性力學、質點力學、剛體力學、約束力學等不同力學條件下的四種體系的MD方法。彈性力學方法是將所考察的原子分子看作剛性球來處理，建立完全彈性碰撞方程，借以求解出原子、分子的運動規律。這種處理可以在液晶的模擬中使用。質點力學模型是將原子、分子作為質點處理，粒子間的相互作用力采用坐標的連續函數。這種力學體系的應用對象非常多，可以用于處理陶瓷、金屬、半導體等無機化合物材料以及有機高分子、生物大分子等幾乎所有的材料。14精選ppt課件2021②力學條件剛體力學方法是把分子作為剛體處理，建立對于剛體的歐拉方程和對于質點系的牛頓方程，聯立求解所研究對象的運動問題，以前主要是在處理像水和四氯化碳那樣的低分子量體系，現已用于研究晶體的相變。約束力學是凍結粒子體系的一部分自由度，進而求解因此而生成約束條件下的質點力學運動方程。對于有機分子來講，因為鍵長、鍵角的振動變化對體系的結構影響不大，將這些自由度凍結是合適的。另外也有采用固定晶體結構的考慮。15精選ppt課件2021③邊界條件問題在處理原子、分子的聚合體問題時，就MD方法而言，能處理的原子(分子)數目要受到計算機運行速度和能力的限制。目前報導的最好水平是能處理109量級的原子數目。這與現實物質含有1023個原子或分子的差距還很大，導致模擬系統原子數少于真實系統的所謂“尺寸效應”的問題。為了減小“尺寸效應”而又不至于使計算工作量過大，對于平衡態MD模擬采用“周期性邊界條件”。16精選ppt課件2021周期性邊界條件圖3-4假設現實物質中的一部分原子(通常為102－105個)，被取出配置在所謂基本單元的箱中，由于基本單元周圍的原子、分子變成表面，從而不同于本來要處理的體狀態。為了防止這種情況，在基本單元周圍配置其復制品(圖3-4)。17精選ppt課件2021周期性邊界條件對于在基本單元周圍邊界的原子、分子所受到的作用力，不僅有來自基本單元之內部的原子或分子的作用力的貢獻，還要考慮來自其近鄰復制單元的原子或分子的作用力的貢獻。這樣，把來自距粒子某一距離（截斷距離）內的j粒子的貢獻（這里與粒子j是否處在基本單元內或復制單元內無關）求和。給出力的方法稱為最小鏡像距離法（MinimumImageDistance）。18精選ppt課件2021周期性邊界條件如果在基本單元內的原子的位置為,周期邊界條件會產生該原子的鏡像。它的位置在

這兒a、b、c是基本單元的三個邊長，l，m，n是整數，取值范圍為-∞到+∞。

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

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

MD模擬給出的直接數據是坐標和速度，對它們進行加工，可得到各種所要的信息。這些信息大致可分為3類：·所模擬系統的平衡態性質·系統處于亞穩態時的結構與性質·系統在遠離平衡態時的動力學過程23精選ppt課件2021平衡態一個物理量A隨時間變化而趨于平衡態的值A0的過程可用下式表示

其中

是馳豫時間。對于小的

，我們只需等待系統達到平衡態后，再開始收集系統在平衡態時的參量。對很大的

，要達平衡態，需要模擬很長的時間，這時MD方法是不合適的。對中等大小的

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

對上式做時間平均，得

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

：,總能量為:

g為總自由度數。標度因子近似可表達成

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

這兒是玻爾茲曼常數，從而得到

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

我們用Einstein關系

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

（3－3－11）Green－Kubo公式和Einstein關系是等價的。32