分子動力學(xué)方法課件

分子動力學(xué)方法分子動力學(xué)方法13.1基本原理計算機模擬分類：(1)隨機模擬方法。優(yōu)點:隨機模擬方法計算的程序簡單，占內(nèi)存少，但是該方法難于處理非平衡態(tài)的問題。(2)分子動力學(xué)方法（MolecularDynamics或簡稱MD）?？梢蕴幚矸瞧胶鈶B(tài)問題。但是使用該方法的程序較復(fù)雜，計算量大，占內(nèi)存也多。分子動力學(xué)方法利用牛頓古典力學(xué)來計算許多分子在相空間中的軌跡。3.1基本原理計算機模擬分類：23.1基本原理分子動力學(xué)(MD)方法的出發(fā)點是物理系統(tǒng)的確定的微觀描述（哈密頓描述方程、拉格朗日方程或牛頓運動方程）。因此，分子動力學(xué)方法是用運動方程來計算多體或少體系統(tǒng)的性質(zhì)，結(jié)果得到的既有系統(tǒng)的靜態(tài)特性，也有動態(tài)特性。蒙特卡羅方法只能得到系統(tǒng)的位形特性。MD方法的具體做法是在計算機上求運動力程的數(shù)值解。為此，通過適當(dāng)?shù)母袷綄Ψ匠踢M行近似，使之適于在計算機上求數(shù)值解。其實質(zhì)是計算一組分子的相空間軌道，其中每個分子各自都服從經(jīng)典運動定律。這里的系統(tǒng)不僅是點粒子系統(tǒng)，也包括具有內(nèi)部結(jié)構(gòu)的粒子組成的系統(tǒng)。3.1基本原理分子動力學(xué)(MD)方法的出發(fā)點是物理系統(tǒng)的確33.1基本原理分子動力學(xué)是在原子、分子水平上求解多體問題的重要的計算機模擬方法。通過求解所有粒子的運動方程，分子動力學(xué)方法可以用于模擬與粒子運動路徑相關(guān)的基本過程。在分子動力學(xué)中，粒子的運動行為是通過經(jīng)典的Newton運動方程所描述。3.1基本原理分子動力學(xué)是在原子、分子水平上求解多體問題的43.1.1粒子運動方程的數(shù)值求解粒子體系的運動方程

Lagrangian方程Lagrangian函數(shù)的定義為

則運動的Lagrangian方程為3.1.1粒子運動方程的數(shù)值求解粒子體系的運動方程Lag53.1.1粒子運動方程的數(shù)值求解粒子的運動方程-Lagrangian方程3.1.1粒子運動方程的數(shù)值求解粒子的運動方程-Lagr63.1.1粒子運動方程的數(shù)值求解單原子的牛頓運動方程3.1.1粒子運動方程的數(shù)值求解單原子的牛頓運動方程73.1.1粒子運動方程的數(shù)值求解Verlet算法

r(t+?t)=r(t)+v(t)?t+(1/2)a(t)?t2(1)

r(t-?t)=r(t)–v(t)?t+(1/2)a(t)?t2(2)

將上面兩式相加，得到：

r(t+?t)=2r(t)–r(t-?t)+a(t)?t2(3)

v(t+?t)=v(t)+a(t)?t+(1/2)b(t)?t2(4)

a(t+?t)=a(t)+b(t)?t(5)

將(5)式的b(t)代入(4)，得到：

v(t+?t)=v(t)+(1/2)[a(t)+a(t+?t)]?t(6)

3.1.1粒子運動方程的數(shù)值求解Verlet算法

83.1.1粒子運動方程的數(shù)值求解其他求解算法：Leap-frogalgorithmr(t+?t)=r(t)+v(t+(1/2)?t)?tv(t+(1/2)?t)=v(t-(1/2)?t)+a(t)?tBeeman’salgorithmr(t+?t)=r(t)+v(t)?t+(2/3)a(t)?t2–(1/6)a(t-?t)?t2

v(t+?t)=v(t)+v(t)?t+(1/3)a(t)?t+(5/6)a(t)?t–

(1/6)a(t?t)?t3.1.1粒子運動方程的數(shù)值求解其他求解算法：93.1.2熱力學(xué)量的計算在物理系統(tǒng)的計算機模擬中，系綜平均必須用時間平均代替，在通常的模擬中，粒子數(shù)N和體積V是固定的。給定初始位置rN(0)和初始動量pN(0)后，一個MD算法將從運動方程生成軌道(rN(t),pN(t))，軌道平均的定義為假定能量守恒，并且軌道在一切具有同一能量的相同休積內(nèi)經(jīng)歷相同的時間，則軌道平均等于微正則系綜平均3.1.2熱力學(xué)量的計算在物理系統(tǒng)的計算機模擬中，系綜平均103.1.2熱力學(xué)量的計算孤立系統(tǒng)的總能量是一個守恒量，沿著分子動力學(xué)模擬生成的任何一條軌道，能量應(yīng)保持不變，即。孤立系統(tǒng)的動能和勢能不是守恒量，它們的大小沿著生成的軌道逐點變化

3.1.2熱力學(xué)量的計算孤立系統(tǒng)的總能量是一個守恒量，沿著113.1.2熱力學(xué)量的計算生成的動能路徑是不連續(xù)的，必須在各個時間間斷點上計算動能的值以求平均:其中從平均動能可以計算系統(tǒng)的溫度，溫度是一個重要的物理量，需要加以監(jiān)測，特別是在模擬的起始階段。3.1.2熱力學(xué)量的計算生成的動能路徑是不連續(xù)的，必須在各123.1.2能量均分定理在熱力學(xué)極限下，一切系綜都是等同的，并且可以應(yīng)用能量均分定理。熱力學(xué)極限是指粒子數(shù)（或體積）趨向無窮大時的極限。一般宏觀物體包含了1023個粒子，可以認(rèn)為是滿足熱力學(xué)極限的。相變只有在熱力學(xué)極限下才會發(fā)生。系統(tǒng)的哈密頓量為則：3.1.2能量均分定理在熱力學(xué)極限下，一切系綜都是等同的，133.1.2能量均分定理由于系統(tǒng)的每個粒子有三個自由度，因此假定位勢在處被截斷。系統(tǒng)內(nèi)部的位形能量的軌道平均值為其中3.1.2能量均分定理由于系統(tǒng)的每個粒子有三個自由度，因此143.1.2能量均分定理由于位勢被截斷，總能量和勢能含有誤差，為了估計必須作出改正。所有的內(nèi)部位形能都加到截止距離為止，尾部改正可以取其中是g(r)對關(guān)聯(lián)函數(shù)，它是在原點r=0處有一個粒子時,在r周圍的體積元dr內(nèi)找到一個粒子的概率，令n(r)為離一個給定粒子的距離在r和r+dr之間的平均粒子數(shù)，于是3.1.2能量均分定理由于位勢被截斷，總能量和勢能含有誤差153.1.2能量均分定理其他的量也需要尾部改正，以壓強的計算為例，這時維里(virial)狀態(tài)方程成立。至于勢能的計算，可以把積分分成兩項，一項是由相互作用力程之內(nèi)的貢獻引起的，一項是對位勢截斷的改正項：長程改正項為3.1.2能量均分定理其他的量也需要尾部改正，以壓強的計算163.2分子動力學(xué)模擬的基本步驟分子動力學(xué)模擬的實際步驟可以劃分為四步：設(shè)定模擬所采用的模型；給定初始條件；趨于平衡的計算過程；宏觀物理量的計算。3.2分子動力學(xué)模擬的基本步驟分子動力學(xué)模擬的實際步驟可以173.2分子動力學(xué)模擬的基本步驟1．模擬模型的設(shè)定硬球勢Lennard-Jones型勢

3.2分子動力學(xué)模擬的基本步驟1．模擬模型的設(shè)定18根據(jù)經(jīng)典物理學(xué)的規(guī)律，可以知道在系綜模擬中的守恒量。微正則系綜的模擬中能量、動量和角動量均為守恒量。在此系綜中他們分別表示為：元胞周期性邊界條件,

3.2分子動力學(xué)模擬的基本步驟根據(jù)經(jīng)典物理學(xué)的規(guī)律，可以知道在系綜模擬中的守恒量。3.219分子動力學(xué)模擬的最小像力約定示意圖

最小像力約定分子動力學(xué)模擬的最小像力約定示意圖最小像力約定202.給定初始條件給定粒子的初始位置和速度的數(shù)值：（1）令初始位置在差分網(wǎng)格格子上，初始速度從玻爾茲曼分布隨機抽樣得到。（2）令初始位置隨機地偏離差分網(wǎng)格格子，初始速度為零。（3）令初始位置隨機地偏離差分網(wǎng)格格子，初始速度從玻爾茲曼分布隨機抽樣得到。3．趨于平衡使系統(tǒng)達到平衡，模擬中需要一個趨衡過程。在這個過程中，我們增加或從系統(tǒng)中移出能量，直到系統(tǒng)具有所要求的能量。3.2分子動力學(xué)模擬的基本步驟2.給定初始條件3.2分子動力學(xué)模擬的基本步驟213.3平衡態(tài)分子動力學(xué)模擬1.微正則系綜的分子動力學(xué)模擬

粒子數(shù)恒定、體積恒定、能量恒定、整個系統(tǒng)的總動量恒等于零。分子動力學(xué)模擬步驟如下（Verlet算法）：（1）給定初始空間位置。（2）計算在第n步時粒子所受的力。（3）計算在第n+1步時所有粒子所處的空間位置。（4）計算第n步的速度。（5）返回到步驟（2），開始下一步的模擬計算。3.3平衡態(tài)分子動力學(xué)模擬1.微正則系綜的分子動力學(xué)模擬22Verlet算法的速度形式：（1）給定初始空間位置。（2）給定初始速度。（3）利用公式：計算在第n+1步時所有粒子所處的空間位置。（4）

計算在第n+1步時所有粒子的速度：（5）返回到步驟（3），開始第n+2步的模擬計算。3.3平衡態(tài)分子動力學(xué)模擬Verlet算法的速度形式：3.3平衡態(tài)分子動力學(xué)模擬233.3平衡態(tài)分子動力學(xué)模擬一般對于能量確定的系統(tǒng)不可能給出精確的初始條件。可以將初始位置固定在格子的格點上，而初始速度則由波爾茲曼分布得出。以上述方法建立的系統(tǒng)一般不會具有所需要的能量，需要在模擬過程中逐漸調(diào)節(jié)系統(tǒng)能量達到給定值。步驟為：(1)對運動方程求解若干步；(2)計算動能和勢能；(3)若能量不等于所需的值，對速度進行標(biāo)度；(4)從第一步開始重復(fù)，直至系統(tǒng)達到平衡為止。3.3平衡態(tài)分子動力學(xué)模擬一般對于能量確定的系統(tǒng)不可能給出24速度標(biāo)度因子:2.正則系綜的分子動力學(xué)模擬速度標(biāo)度因子:

3.3平衡態(tài)分子動力學(xué)模擬速度標(biāo)度因子:3.3平衡態(tài)分子動力學(xué)模擬253.3平衡態(tài)分子動力學(xué)模擬正則系綜分子動力學(xué)的模擬具體步驟:(Verlet算法的速度形式)3.3平衡態(tài)分子動力學(xué)模擬正則系綜分子動力學(xué)的模擬26（1）給定初始空間位置，（2）給定初始速度，（3）利用公式：計算在第n+1步時所有粒子所處的空間位置，（4）計算在第n+1步時所有粒子的速度：

動能和速度標(biāo)度因子：

（5）計算將速度乘以標(biāo)度因子的值，并讓該值作為下一次計算時，第n+1步粒子的速度：

（6）

返回到步驟（3），開始第n+2步的模擬計算。

（1）給定初始空間位置，273.4MD在材料科學(xué)中的應(yīng)用3.4.1高分子鏈動力學(xué)模擬計明娟等用高溫淬火分子動力學(xué)模擬方法研究了甲硫氨酸一腦啡肽(Met-enkephalin)在真空中的構(gòu)象性質(zhì)，經(jīng)聚類分析和能量優(yōu)化得到了十三個低能構(gòu)象，并與嗎啡進行了空間擬合。3.4MD在材料科學(xué)中的應(yīng)用3.4.1高分子鏈動力學(xué)模擬283.4MD在材料科學(xué)中的應(yīng)用按圖l的序列構(gòu)造出Met-enkephalin分子并進行初始的優(yōu)化，然后以每0.1ps溫度增加l0K進行動力學(xué)運算，加熱到800K。在此溫度下進行高溫淬火分子動力學(xué)模擬計算。平衡時間為10ps，模擬時間為30ps，時間步長為1.0fs，取樣間隔為0.2ps，在正則(NTV)系綜下運行．對平衡后的構(gòu)象作進一步分析，根據(jù)它們骨架(φ，Ψ)構(gòu)象分布的分析，以及所有結(jié)構(gòu)中每一個殘基Co-C骨架片段之間均方根偏差不同性的計算產(chǎn)生五個等級，將均方根偏差最小間隔(0.1nm≤rms≤0.2nm)的結(jié)構(gòu)分成十三個簇，從每個簇中選取最低能量構(gòu)象作為此簇代表，最后將十三個所選結(jié)構(gòu)中的每一個結(jié)構(gòu)利用最陡下降法和共軛梯度法分別運行1000步進行完全的能量優(yōu)化。直至rms能量偏差<0.004kJ?mol而收斂，然后將優(yōu)化后的十三個構(gòu)象與嗎啡進行空間擬合。3.4MD在材料科學(xué)中的應(yīng)用按圖l的序列構(gòu)造出Met-en29圖3.2經(jīng)QMD模擬所得的十三個低能構(gòu)象圖3.2經(jīng)QMD模擬所得的十三個低能構(gòu)象303.4MD在材料科學(xué)中的應(yīng)用模擬結(jié)果表明，Met-enkephalin是一種構(gòu)象柔性大的多肽分子，可具有多種構(gòu)象，主要為Gly2-Gly3和Gly3-Phe4之間的β折疊型式。大多數(shù)的結(jié)構(gòu)具有由氫鍵連接兩個末端殘基所形成的假大環(huán)，側(cè)鏈和芳香環(huán)定位在假大環(huán)的同側(cè)。酪氨酸殘基中的氨基N相當(dāng)于嗎啡中哌啶環(huán)上的季胺N，酪氨酸殘基中的酚環(huán)相當(dāng)于嗎啡中的酚環(huán)，苯丙氨酸殘基中的芳香環(huán)相當(dāng)于嗎啡中環(huán)己烯環(huán)，此肽與嗎啡在空間結(jié)構(gòu)(藥效基團)上非常相似，都滿足阿片受體的要求。所以，它們可作用于同一受體。3.4MD在材料科學(xué)中的應(yīng)用模擬結(jié)果表明，Met-enke313.4MD在材料科學(xué)中的應(yīng)用3.4.2脆性斷裂模擬單德彬等為了研究在微觀尺度上單晶銅的彎曲裂紋萌生和擴展機理，建立了單晶銅彎曲變形的分子動力學(xué)模型，應(yīng)用速度標(biāo)定法控制溫度，采用Morse勢進行了單晶銅的彎曲變形分子動力學(xué)模擬。3.4MD在材料科學(xué)中的應(yīng)用3.4.2脆性斷裂模擬323.4MD在材料科學(xué)中的應(yīng)用二維分子動力學(xué)模型如圖3.3。圖3.3彎曲變形仿真實驗?zāi)Ｐ?無凹槽模型)圖3.4彎曲變形仿真實驗?zāi)Ｐ?無凹槽模型)3.4MD在材料科學(xué)中的應(yīng)用二維分子動力學(xué)模型如圖3.3。333.4MD在材料科學(xué)中的應(yīng)用以上模擬晶面是(100)晶面，原子位置按理想點陣排列。原子運動按牛頓原子處理，牛頓原子的作用由Morse勢函數(shù)計算。仿真實施后分別在模型的左右兩端3列原子上施加向下的位移，同時施加彎矩(M)。計算采用速度-Verlet算法。原子的初始位置設(shè)在fcc晶體的晶格上，原子的初始速度可利用Maxwell分布確定，也可直接置為零。溫度控制采用速度標(biāo)定法控制在絕對零度，以避免原子熱激活的復(fù)雜影響。3.4MD在材料科學(xué)中的應(yīng)用以上模擬晶面是(100)晶面，34圖3.5單晶Cu(100)晶面無凹槽模型彎曲變形仿真實驗圖3.5單晶Cu(100)晶面無凹槽模型彎曲變形仿真實驗35圖3.6單晶Cu(100)晶面有凹槽模型彎曲變形仿真實驗圖3.6單晶Cu(100)晶面有凹槽模型彎曲變形仿真實驗363.4MD在材料科學(xué)中的應(yīng)用通過模擬試驗結(jié)果圖3.5發(fā)現(xiàn)，在位移和彎矩的作用下，個別原子間間隙的不斷增大和應(yīng)變能的不斷積累最終使晶體內(nèi)部出現(xiàn)空位，在128ps已經(jīng)很明顯的看到空位，這些空位位于模型中部，空位的不斷合并萌生裂紋，微觀裂紋的頂端又形成空位，隨著時間步的推移，塑性變形加劇，裂紋逐漸萌生和擴展，發(fā)展成為貫穿模型的裂紋，形成裂紋的擴展，微觀裂紋在裂紋頂端形核而造成主裂紋的連續(xù)擴展，裂尖向下萌生，后續(xù)微觀裂紋的擴展類似于宏觀裂紋擴展，與無凹槽模型相比，圖3.6有凹槽模型其裂紋萌生和擴展的時間明顯縮短，表明裂紋缺陷對斷裂過程起促進作用。3.4MD在材料科學(xué)中的應(yīng)用通過模擬試驗結(jié)果圖3.5發(fā)現(xiàn)，373.4MD在材料科學(xué)中的應(yīng)用研究結(jié)果表明：應(yīng)變能的不斷積累使晶體內(nèi)部產(chǎn)生空位，材料的裂紋萌生于空位，空位的合并形成納米級裂紋，后續(xù)微觀裂紋的擴展類似于宏觀裂紋；裂紋缺陷促進了裂紋的萌生和擴展。而且對無凹槽模型和有凹槽模型裂紋萌生和擴展的時間步長的對比研究表明，裂紋缺陷對斷裂過程起促進作用。3.4MD在材料科學(xué)中的應(yīng)用研究結(jié)果表明：應(yīng)變能的不斷積累383.4MD在材料科學(xué)中的應(yīng)用鐵中裂紋擴展的結(jié)構(gòu)演化吳映飛用分子動力學(xué)模擬的方法研究了bcc—Fe中I型裂紋在應(yīng)力及溫度場下裂尖區(qū)的結(jié)構(gòu)演化問題。圖3.7裂紋的幾何結(jié)構(gòu)示意圖

3.4MD在材料科學(xué)中的應(yīng)用鐵中裂紋擴展的結(jié)構(gòu)演化圖3.7393.4MD在材料科學(xué)中的應(yīng)用裂紋前沿沿[101]方向，裂紋面垂直于[0l0]方向，裂紋擴展沿[101]。加載通過增加應(yīng)力強度因子(KⅠ)實現(xiàn)，K1由臨界應(yīng)力強度因子KIC度量?；诟飨蛲跃€彈性力學(xué)的平面應(yīng)變條件得到初始裂紋，進行了兩組模擬。(1)第一組模擬選取初始溫度T=5K、應(yīng)力強度因子KⅠ=1.0KIC，設(shè)在所對應(yīng)的外載下存在初始裂紋，隨后溫度增加到l00，300，和500K，外載增加到KⅠ=2.8KIC。在這組模擬中，顯示出堆垛層錯，裂紋形狀尖銳，且擴展速度較快，(2)第二組模擬選取初始溫度T=100K，應(yīng)力場條件同于第一組，溫度加到300和500K，外載加到KⅠ=2.8KIC。在這組模擬中，顯示出了堆垛層錯、位錯發(fā)射、裂尖鈍化和分枝，以及孿晶帶；該組裂紋的擴展速度慢于第一組。

3.4MD在材料科學(xué)中的應(yīng)用裂紋前沿沿[101]方向，裂紋40圖3.9KI=1.8KIC后，裂尖區(qū)出現(xiàn)了孿晶帶(T=100K,KI=2.5KIC)圖3.8低溫下裂尖區(qū)原子的ABAB型堆垛層錯(T=5K,KI=1.8KIC)圖3.9KI=1.8KIC后，裂尖區(qū)出現(xiàn)了孿晶帶圖3.8413.4MD在材料科學(xué)中的應(yīng)用研究結(jié)果表明：在裂紋擴展過程中，裂尖區(qū)結(jié)構(gòu)演化與初始裂紋的工作條件相關(guān)，并且依賴于溫度與外載的協(xié)同作用。在相同條件下裂紋周圍局域原子的能量也受初始裂紋的工作條件影響。3.4MD在材料科學(xué)中的應(yīng)用研究結(jié)果表明：在裂紋擴展過程中423.4MD在材料科學(xué)中的應(yīng)用3.4.3表面效應(yīng)對于晶界附近原子遷移的影響K.Saitoh等采用分子動力學(xué)模擬方法（對微電子器件中鋁芯絲材內(nèi)晶界附近原子的遷移進行了分析。該模擬是基于EMT理論建立的模型，這一理論結(jié)合大量電子濃度的體效應(yīng)可以驗證出更加準(zhǔn)確的晶界或表面處原子結(jié)構(gòu)的研究結(jié)果。主要工作分為兩部分，即晶界附近擴散性能的模擬和晶界附近原子的重新排列。它的算法是對MD時間段（2.0fs）和運動方程的積分，在該模型中，將∑=5（100）晶界和自由表面分別設(shè)定為xz面和yz面。3.4MD在材料科學(xué)中的應(yīng)用3.4.3表面效應(yīng)對于晶界附近43圖3.10擴散性能研究中的原子模型圖3.10擴散性能研究中的原子