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

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

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

4、見表3.1。年代創立者創造內容工作（MD分類名稱）1957B. J Alder & T. E. Wainwright剛性球MD方法1963A. Rahman質點系MD方法1971Rahman &F. H. Stillinger剛性系統MD方法1972W. Lees &S. F. Edwards平衡系統MD方法（存在速度梯度）1977J. P. Rychaert et al.約束系統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巨正則系統MD方法表3.1 MD方法的里程碑工作方法的里程碑工作 3.2 MD方法計算初步方法計算初步 在計算機出現以前，作為根據原子間相互作用力等微觀信息了解多原子或分子團的結構、性質的方法，所采用的是基于統計理論的數學解析法。然而，原子間相互作用力稍微復雜一些，不用說求解統計理論嚴格方程解，就是進行數值求解也是一件很困難的事。 MD方法就是數值求解多體系統的確定性運動方程，并

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

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

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

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

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

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

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

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

14、GreenKubo formula 擴散系數D也可用速度自相關函數通過下面的GreenKubo 公式表示 （3311） GreenKubo公式和Einstein 關系是等價的。 101133NiiiDv tv odtvtvo dtNPartial diffusion 在多種原子組成的體系中（例如AB2），可以對某一種原子求和，得到Partial diffusion，從而得到原子A或原子B的擴散系數DA、DB。非均勻體系中的擴散非均勻體系中的擴散 對一非均勻體系，例如一塊材料，有體和晶界部分。可分別對某一類原子做平均，如對晶界原子，或對體中的原子，得到不同部分原子不同的擴散。與溫度有關的擴散與溫

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