第6章 分子動(dòng)力學(xué)方法

1、- 28 - 第6章 分子動(dòng)力學(xué)方法 第6章 分子動(dòng)力學(xué)方法經(jīng)典分子動(dòng)力學(xué)方法無疑是材料，尤其是大分子體系和大體系模擬有效的方法之一。分子動(dòng)力學(xué)可以用于NPT，NVE，NVT等不同系綜的計(jì)算，是一種基于牛頓力學(xué)確定論的熱力學(xué)計(jì)算方法。與蒙特卡羅法相比在宏觀性質(zhì)計(jì)算上具有更高的準(zhǔn)確度和有效性，可以廣泛應(yīng)用于物理，化學(xué)，生物，材料，醫(yī)學(xué)等各個(gè)領(lǐng)域。本章在介紹分子動(dòng)力學(xué)的基本概念的基礎(chǔ)上，簡單介紹了分子動(dòng)力學(xué)的基本思想，勢函數(shù)分類和基本方程。然后介紹了分子動(dòng)力學(xué)的常用系綜和典型的NPT，NVE，NVT系綜基本方程。結(jié)合材料建模中的基本簡化方法和技巧，闡述了邊界條件和時(shí)間積分的數(shù)值處理技巧。最后，利用

2、統(tǒng)計(jì)力學(xué)的基本概念給出分子動(dòng)力學(xué)的計(jì)算信息的解析方式。并且結(jié)合Materials Explore軟件計(jì)算分析了CNT的幾何結(jié)構(gòu)穩(wěn)定性。6.1 引言分子動(dòng)力學(xué)方法(Molecular Dynamics, MD)方法是一種按該體系內(nèi)部的內(nèi)稟動(dòng)力學(xué)規(guī)律來計(jì)算并確定位形的變化的確定性模擬方法。首先需要在給定的外界條件下建立一組粒子的運(yùn)動(dòng)方程，然后通過直接對系統(tǒng)中的一個(gè)個(gè)粒子運(yùn)動(dòng)方程進(jìn)行數(shù)值求解，得到每個(gè)時(shí)刻各個(gè)分子的坐標(biāo)與動(dòng)量，即在相空間的運(yùn)動(dòng)軌跡，再利用統(tǒng)計(jì)力學(xué)方法得到多體系統(tǒng)的靜態(tài)和動(dòng)態(tài)特性，從而獲得系統(tǒng)的宏觀性質(zhì)?？梢钥闯?，分子動(dòng)力學(xué)方法中不存在任何隨機(jī)因素，這個(gè)也是分子動(dòng)力學(xué)方法和后文要提到的

3、蒙特卡洛方法的區(qū)別之一。在分子動(dòng)力學(xué)方法的處理過程中，方程組的建立是通過對物理體系的微觀數(shù)學(xué)描述給出的。在這個(gè)微觀的物理體系中，每個(gè)分子都各自服從經(jīng)典的牛頓力學(xué)定律(或者是拉格朗日方程)。每個(gè)分子運(yùn)動(dòng)的內(nèi)稟動(dòng)力學(xué)是用理論力學(xué)上的哈密頓量或者拉格朗日函數(shù)來描述，也可以直接用牛頓運(yùn)動(dòng)方程來描述。確定性方法是實(shí)現(xiàn)玻爾茲曼的統(tǒng)計(jì)力學(xué)途徑。這種方法可以處理與時(shí)間有關(guān)的過程，因而可以處理非平衡態(tài)問題。但是分子動(dòng)力學(xué)方法的計(jì)算機(jī)程序相對蒙特卡羅較復(fù)雜，其計(jì)算成本較高。分子動(dòng)力學(xué)方法發(fā)展歷史改革經(jīng)歷了近60年，分子動(dòng)力學(xué)方法是20世紀(jì)50年代后期由Alder B J 和Wainwright T E創(chuàng)造發(fā)展的。

4、Alder和Wainwright在1957年利用分子動(dòng)力學(xué)模擬，發(fā)現(xiàn)了“剛性球組成的集合系統(tǒng)會(huì)發(fā)生由其液相到結(jié)晶相的相轉(zhuǎn)變”，后來人們稱這種相變?yōu)锳lder相變。其結(jié)果表明，不具有引力的粒子系統(tǒng)也具有凝聚態(tài)。到20世紀(jì)70年代，產(chǎn)生了剛性體系的動(dòng)力學(xué)方法，成功地被應(yīng)用于水和氮等分子性溶液體系的處理；1972年，Less A W和 Edwards S F 等發(fā)展了該方法并擴(kuò)展到了存在速度梯度，即處于非平衡狀態(tài)的系統(tǒng)。之后，此方法被 Ginan M J等推廣到了具有溫度梯度的非平衡系統(tǒng)，從而構(gòu)造并形成了所謂的非平衡分子動(dòng)力學(xué)方法體系。進(jìn)人20世紀(jì)80年代之后，出現(xiàn)了在分子內(nèi)部對一部分自由度施加約束

5、條件的新的分子動(dòng)力學(xué)方法，從而使 分子動(dòng)力學(xué)方法可適用于類似蛋白質(zhì)等生物大分子的解析與設(shè)計(jì)。分子動(dòng)力學(xué)方法真正作為材料科學(xué)領(lǐng)域的一個(gè)重要研究方法，開始于由Andersen, Parrinello和Rahman等創(chuàng)立恒壓分子動(dòng)力學(xué)方法和Nse等完成恒溫分子動(dòng)力學(xué)方法的建立及在應(yīng)用方面的成功。后來，針對勢函數(shù)模型化比較困難的半導(dǎo)體和金屬等，1985 年人們又提出了將基于密度泛函理論的電子論和分子動(dòng)力學(xué)方法有機(jī)統(tǒng)一起來的所謂Car-Parrinello 方法，亦即第一性原理分子動(dòng)力學(xué)方法。這樣，分子動(dòng)力學(xué)的方法進(jìn)一步得到發(fā)展和完善，它不僅可以處理半導(dǎo)體和金屬的間題，同時(shí)還可應(yīng)用于處理有機(jī)物和化學(xué)反應(yīng)

6、。關(guān)于Car-Parrinello分子動(dòng)力學(xué)方法將在第7章重點(diǎn)學(xué)習(xí)討論。本章重點(diǎn)討論經(jīng)典分子動(dòng)力學(xué)方法的基本原理和計(jì)算方法。6.2 分子動(dòng)力學(xué)的計(jì)算框架6.2.1 基本思想分子動(dòng)力學(xué)方法是沿用牛頓運(yùn)動(dòng)方程或拉格朗日方程、哈密頓方程，通過考察粒子的運(yùn)動(dòng)來研究多粒子系統(tǒng)的物理性質(zhì)。在處理孤立粒子或原子團(tuán)簇(不考慮與其他粒子的相互作用)時(shí)，單純地對牛頓運(yùn)動(dòng)方程進(jìn)行積分求解就行了；而在處理凝聚系統(tǒng)時(shí)，可以認(rèn)為將所考慮的N個(gè)粒子放入具有一定體積的容器中，在這樣組成的封閉系中，建立直角坐標(biāo)(x, y, z)，每個(gè)粒子的位置就由三個(gè)坐標(biāo)分量決定，通過求解3N個(gè)聯(lián)立方程組就可以得到保守系的總能量。分子動(dòng)力學(xué)的

7、基本思想在于通過設(shè)定原子之間的相互作用(勢函數(shù))和相關(guān)的系綜(亦即作用對象和條件)，確定其基本的模擬范疇。在給定的牛頓方程、拉格朗日方程或者是哈密頓方程進(jìn)行時(shí)間的迭代，在達(dá)到指定的力學(xué)收斂條件之后得到一個(gè)最終的位置坐標(biāo)，然后通過相關(guān)的位形信息和運(yùn)動(dòng)的速度和加速度等信息通過熱力統(tǒng)計(jì)學(xué)的方法，給出模擬體系的統(tǒng)計(jì)信息，主要包括復(fù)雜的光學(xué)，電學(xué)和一些其他的物理信息。其基本思想，總結(jié)如圖6-1所示。圖 6-1 分子動(dòng)力學(xué)的一般思想6.2.2 計(jì)算流程總體來說，分子動(dòng)力學(xué)的基本計(jì)算有如下四個(gè)步驟：1. 模型的設(shè)定，也就是勢函數(shù)的選取勢函數(shù)的研究和物理系統(tǒng)中對物質(zhì)的描述研究息息相關(guān)。分子動(dòng)力學(xué)的模擬最早使用

8、是硬球勢，即小于臨界值時(shí)無窮大，大于等于臨界值時(shí)為零。常用的是Lennard-Jones勢函數(shù)，還有嵌入原子EAM勢函數(shù)，不同的物質(zhì)狀態(tài)描述用不同的勢函數(shù)。模型勢函數(shù)一旦確定，就可以根據(jù)物理學(xué)規(guī)律求得模擬中的守恒量。關(guān)于勢函數(shù)選取將在6.4節(jié)進(jìn)一步討論。2. 給定初始條件運(yùn)動(dòng)方程的求解需要知道粒子的初始位置和速度，不同的算法要求不同的初始條件。如Verlet算法需要兩組坐標(biāo)來啟動(dòng)計(jì)算，一組零時(shí)刻的坐標(biāo)，一組是前進(jìn)一個(gè)時(shí)間步的坐標(biāo)或者一組零時(shí)刻的速度值。一般意義上講系統(tǒng)的初始條件不可能知道，實(shí)際上也不需要精確選擇代求系統(tǒng)的初始條件，因?yàn)楫?dāng)模擬時(shí)間足夠長時(shí)，系統(tǒng)就會(huì)忘掉初始條件(對于無記憶的體系而

9、言)。當(dāng)然，合理的初始條件可以加快系統(tǒng)趨于平衡的時(shí)間和步程，獲得好的精度。常用的初始條件可以選擇為，令初始位置在差分劃分網(wǎng)格的格子上，初始速度則從玻爾茲曼分布隨機(jī)抽樣得到；令初始位置隨機(jī)的偏離差分劃分網(wǎng)格的格子上，初始速度為零；令初始位置隨機(jī)的偏離差分劃分網(wǎng)格的格子上，初始速度也是從玻爾茲曼分布隨機(jī)抽樣得到。圖6-2 分子動(dòng)力學(xué)計(jì)算流程3. 趨于平衡計(jì)算在邊界條件和初始條件給定后就可以解運(yùn)動(dòng)方程，進(jìn)行分子動(dòng)力學(xué)模擬。但這樣計(jì)算出的系統(tǒng)是不會(huì)具有所要求的系統(tǒng)能量，并且這個(gè)狀態(tài)本身也還不是一個(gè)平衡態(tài)。為使得系統(tǒng)平衡，模擬中設(shè)計(jì)一個(gè)趨衡過程，即在這個(gè)過程中，增加或者從系統(tǒng)中移出能量，直到持續(xù)給出確定

10、的能量值，稱此時(shí)系統(tǒng)已經(jīng)達(dá)到平衡，達(dá)到平衡的時(shí)間稱為弛豫時(shí)間。 分子動(dòng)力學(xué)中，積分步程的大小選擇十分重要，這決定了模擬所需要的時(shí)間。為了減小誤差，步長要小，但過小系統(tǒng)模擬的弛豫時(shí)間就長。因此根據(jù)經(jīng)驗(yàn)選擇適當(dāng)?shù)牟介L。如，對一個(gè)具有幾百個(gè)氬氣分子的體系，采用Lennard-Jones勢函數(shù)，發(fā)現(xiàn)取h為0.01量級，可以得到很好的結(jié)果。這里選擇的h是沒有量綱的，實(shí)際上這樣選擇的h對應(yīng)的時(shí)間在10-14s的量級。如果模擬1000步，系統(tǒng)達(dá)到平衡，弛豫時(shí)間只有10-11s。4. 宏觀物理量的計(jì)算實(shí)際計(jì)算宏觀的物理量往往是在模擬的最后階段進(jìn)行的，它是沿相空間軌跡求平均來計(jì)算得到的，根據(jù)各態(tài)歷經(jīng)假說，時(shí)間平

11、均代替系綜平均。關(guān)于各態(tài)歷經(jīng)假說將在第8章的蒙特卡洛方法中深入討論。6.2.3 經(jīng)典分子動(dòng)力學(xué)中近似處理實(shí)際的計(jì)算體系由于其外在條件的復(fù)雜性，往往需要對計(jì)算的系統(tǒng)進(jìn)行一定的近似處理，在經(jīng)典分子動(dòng)力學(xué)計(jì)算中，引進(jìn)以下近似處理：v 經(jīng)典粒子相互作用，不考慮電子相互作用量子效應(yīng)。v 力的作用形式，由參數(shù)可調(diào)的相互作用勢函數(shù)決定，并經(jīng)過實(shí)驗(yàn)測定來驗(yàn)證。v 模擬體系與實(shí)際體系相差較大，一般需要采用周期邊界來擴(kuò)展計(jì)算體系。時(shí)間平均是在有限時(shí)間內(nèi)完成?！揪毩?xí)與思考】6-1. 查找文獻(xiàn)，根據(jù)上述的分子動(dòng)力學(xué)的處理流程圖，編寫實(shí)現(xiàn)分子動(dòng)力學(xué)的簡單程序，可參考Daan F和Berend S編著的分子模擬一書。6.

12、3 分子動(dòng)力學(xué)的系綜系綜(ensemble)是指在一定的宏觀條件下(約束條件)，大量性質(zhì)和結(jié)構(gòu)完全相同的、處于各種運(yùn)動(dòng)狀態(tài)的、各自獨(dú)立的系統(tǒng)集合，或稱為統(tǒng)計(jì)系綜。系綜是用統(tǒng)計(jì)方法描述熱力學(xué)系統(tǒng)的統(tǒng)計(jì)規(guī)律性時(shí)引入的一個(gè)基本概念；系綜是統(tǒng)計(jì)理論的一種表述方式，系綜理論使統(tǒng)計(jì)物理成為普遍的微觀統(tǒng)計(jì)理論；系綜并不是實(shí)際的物體，構(gòu)成系綜的系統(tǒng)才是實(shí)際物體。 約束條件是由一組外加宏觀參量來表示。在平衡統(tǒng)計(jì)力學(xué)范疇下，可以用來處理穩(wěn)定系綜。6.3.1 常用系綜分類 1. 正則系綜(canonical ensemble)全稱應(yīng)為“宏觀正則系綜”，簡寫為NVT，即表示具有確定的粒子數(shù)(N)、體積(V)、溫度(T

13、)。正則系綜是分子動(dòng)力學(xué)方法模擬處理的典型代表。假定N個(gè)粒子處在體積為V的盒子內(nèi)，將其埋入溫度恒為T的熱浴中。此時(shí)，總能量(E)和系統(tǒng)壓強(qiáng)(P)可能在某一平均值附近起伏變化。平衡體系代表封閉系統(tǒng)，與大熱源熱接觸平衡的恒溫系統(tǒng)。正則系綜的特征函數(shù)是亥姆霍茲自由能。 2. 微正則系綜(micro-canonical ensemble)簡寫為NVE，即表示具有確定的粒子數(shù)(N)、體積(V)、總能量(E)。微正則系綜廣泛被應(yīng)用在分子動(dòng)力學(xué)模擬中。假定N個(gè)粒子處在體積為V的盒子內(nèi)，并固定總能量(E)。此時(shí)，系綜的溫度(T)和系統(tǒng)壓強(qiáng)(P)可能在某一平均值附近起伏變化。平衡體系為孤立系統(tǒng)，與外界即無能量交

14、換，也無粒子交換。微正則系綜的特征函數(shù)是熵。3. 等溫等壓(constant-pressure, constant-temperature)簡寫為NPT，即表示具有確定的粒子數(shù)(N)、壓強(qiáng)(P)、溫度(T)。其總能量(E)和系統(tǒng)體積(V)可能存在起伏。體系是可移動(dòng)系統(tǒng)壁情況下的恒溫?zé)嵩?。特征函?shù)是吉布斯自由能。 4. 等壓等焓(constant-pressure, constant-enthalpy)簡寫為NPH，即表示具有確定的粒子數(shù)(N)、壓強(qiáng)(P)、焓(H)。由于，故在該系綜下進(jìn)行模擬時(shí)要保持壓力與焓值為固定，其調(diào)節(jié)技術(shù)的實(shí)現(xiàn)也有一定的難度，這種系綜在實(shí)際的分子動(dòng)力學(xué)模擬中已經(jīng)很少遇到了

15、。 5. 巨正則系綜(grand canonical ensemble)簡寫為VT，即表示具有確定的粒體積(V)、溫度(T)和化學(xué)勢()。巨正則系綜通常是蒙特卡羅模擬的對象和手段。此時(shí)、系統(tǒng)能量(E)、壓強(qiáng)(P)和粒子數(shù)(N)會(huì)在某一平均值附近有一個(gè)起伏。體系是一個(gè)開放系統(tǒng)，與大熱源大粒子源熱接觸平衡而具有恒定的溫度(T)。特征函數(shù)是馬休(Massieu)函數(shù)J(，V，T)。對于不同的系綜，由于其選擇的運(yùn)動(dòng)方程的求解方法亦有不同。針對不同的計(jì)算體系具體的細(xì)節(jié)比較復(fù)雜，下文將以NEV和NTP這兩類簡單的系綜運(yùn)用兩類不同的條件方法下求解的不同方程做簡要的介紹。6.3.2 NEV系綜基本方程N(yùn)EV系

16、綜是一類較為常見的系綜。對于該系綜，經(jīng)典分子動(dòng)力學(xué)方法一般采用差分近似法求解牛頓運(yùn)動(dòng)方程，并追蹤系統(tǒng)的時(shí)間變化?？紤]由N個(gè)粒子構(gòu)成的系統(tǒng)，假設(shè)第i個(gè)粒子在時(shí)刻位置，速度為受到的作用力為，為求出該粒子(i)在時(shí)刻的位置，經(jīng)常使用的方法是Verlet算法(對于該算法的具體細(xì)節(jié)在下文將有詳細(xì)的推導(dǎo))。根據(jù)系統(tǒng)的動(dòng)能及統(tǒng)計(jì)物理學(xué)中的能量均分定理，則系統(tǒng)的平衡溫度可表達(dá)為(6.1)根據(jù)此方法，在晶體中粒子的配置與排布、晶形變化等結(jié)構(gòu)相變也可以用計(jì)算機(jī)模擬進(jìn)行研究。根據(jù)物理化學(xué)的基本氣體的基本方程可以知道，體系的NEV系綜下的計(jì)算體系主要變化的量有兩個(gè)壓力和溫度。下文就從恒溫和恒壓兩個(gè)條件下對NEV系綜的

17、基本方法做簡單的推導(dǎo)。1. 恒溫方法Nóse S和Hoover W G引入與熱源相關(guān)的參數(shù)，嘗試構(gòu)造恒溫狀態(tài)，亦即具有確定N，V，T值的系統(tǒng)，對此可設(shè)想為與大熱源接觸而達(dá)到平衡的系統(tǒng)。由于熱源很大，交換能量不會(huì)改變熱源的溫度。在兩者建立平衡后，系統(tǒng)將與熱源具有相同的溫度，系統(tǒng)與熱源合起來構(gòu)成一個(gè)復(fù)合系統(tǒng)，如圖6-3所示。這樣，系統(tǒng)中微觀粒子的動(dòng)力學(xué)方程：(6.2)(6.3)(6.4)式中為系統(tǒng)的勢函數(shù)，為有效質(zhì)量，表示熱力學(xué)摩爾系數(shù)。式(6.2)至式(6.4)同往常的分子動(dòng)力學(xué)方法的區(qū)別體現(xiàn)在式(6.3)中，即增加了與熱源的相互作用相關(guān)的并與力的量綱相同的一項(xiàng)()。與熱源相關(guān)的變化參

18、數(shù)的運(yùn)動(dòng)方程表明，當(dāng)系統(tǒng)的總能量大于時(shí)，是增加的，從而顯示出使粒子速度減小那樣的作用；反之則顯示出使粒子速度增大；是表示與溫度控制有關(guān)的常數(shù)。圖6-3 采取擴(kuò)展系方法恒溫分子動(dòng)力學(xué)方法原理示意圖該方法的特征是系統(tǒng)處于微觀狀態(tài)的分布服從正則分布，若在推導(dǎo)式(6.2)、式(6.3)、式(6.4)中，引人由公式定義的變量，則系統(tǒng)的哈密頓量(即能量)可以表示為(6.5)式中第二項(xiàng)和第三項(xiàng)對應(yīng)于與熱源相聯(lián)系的部分。這說明，由于系統(tǒng)與熱源存在熱接觸，二者可以交換能量，因此粒子的能量發(fā)生變化(H0H)，同時(shí)也說明系統(tǒng)可能的微觀狀態(tài)可具有不同的能量值。2. 恒壓方法為了進(jìn)行壓力的控制，考慮如圖6-4所示的方法

19、，即利用活塞調(diào)控系統(tǒng)的體積變化。對于質(zhì)量一定的理想氣體系統(tǒng)，體積變大，其壓力就變小，反之亦然。利用這樣的原理，可以通過改變體積使壓力達(dá)到某一個(gè)預(yù)定值。就實(shí)際的模擬而言，可控制體積大小均勻(三維)地變化，將粒子的坐標(biāo)、動(dòng)量表示成如下形式(6.6)從而將粒子的坐標(biāo)、動(dòng)量與對分子動(dòng)力學(xué)單元(邊長)歸一化的參變量和聯(lián)系起來。應(yīng)用歸一化的參變量并根據(jù)正則方程可得到粒子運(yùn)動(dòng)方程。圖6-4 恒壓分子動(dòng)力學(xué)方法原理示意圖(實(shí)際上，活塞為三維運(yùn)動(dòng))粒子系和活塞組成復(fù)合系統(tǒng)，其哈密頓量可由下式給出(6.7)(6.8)式中，是外壓；是體積對應(yīng)的共軛動(dòng)量；是決定于體積變化速度的因子。若將由上述哈密頓量導(dǎo)出的運(yùn)動(dòng)方程用

20、粒子的實(shí)際坐標(biāo)直接寫出來，則有(6.9)(6.10)(6.11)式(6.9)中右邊第一項(xiàng)是根據(jù)維里(Virial)定理求出的系統(tǒng)的瞬間壓P。在此處理中，特別重要的是式(6.9)右邊第一項(xiàng)與統(tǒng)計(jì)力學(xué)壓力表達(dá)式相同。6.3.3 NTP 系綜質(zhì)點(diǎn)系的基本方程前面已討論了關(guān)于NEV系綜的質(zhì)點(diǎn)系方程。鑒于NTV和NHP系綜的基本方程包含于NTP系綜的基本方程之中，在此只對NTP系綜的基本方程作一些說明。NTP系綜是在宏觀上具有確定的粒子數(shù)N、溫度T和壓力P的系綜。在此系綜中，T和P可作為外部參數(shù)給出。那么，怎樣才能滿足N，T，P恒定的條件呢？顯然，粒子數(shù)N保持一定是最簡單的，因?yàn)橹灰顾疾斓幕締卧獌?nèi)

21、的粒子數(shù)恒定就行了。1. 壓力恒定體系為簡單起見，以由圓柱容器內(nèi)充滿氣體和調(diào)節(jié)活塞構(gòu)成的系統(tǒng)為例，如圖6-5所示。設(shè)活塞的重量為W，則在平衡狀態(tài)下，由活塞的重量產(chǎn)生的外部壓力()與由氣體產(chǎn)生的內(nèi)壓平衡。如果外壓變大，為了保持平衡，則氣體部分的體積將收縮，從而內(nèi)壓變大。圖6-5 活塞與氣體容器構(gòu)成的示意圖因此，把體積作為系統(tǒng)的動(dòng)力學(xué)參數(shù)建立其運(yùn)動(dòng)方程，就好像維持壓力恒定。其方程可寫為(6.12)式中，為體積，是由維里(Virial)定理求出的內(nèi)壓，是對應(yīng)于的假想質(zhì)量。圖6-5為活塞與氣體容器構(gòu)成的示意圖。將上述想法推廣到一般情況下的晶格系統(tǒng)，如圖6-6所示，對一般的晶體，其基本單元可由三個(gè)平移矢

22、量 , 決定其形狀和大小。將的三個(gè)分量寫成矩陣，該矩陣規(guī)定了晶體基本單元的形狀和大小，這樣也成為原子的實(shí)際坐標(biāo)和晶格矢量之間的變換矩陣，亦即圖6-6 基于單元形狀示意圖在式(6.15)中，體積作為動(dòng)力學(xué)參數(shù)，而在一般晶體系統(tǒng)中，應(yīng)該作為動(dòng)力學(xué)參數(shù)。若給出此動(dòng)力學(xué)系統(tǒng)的哈密頓量，則有(6.13)式中，是與格矢相應(yīng)的原子共扼動(dòng)量，為數(shù)值行列式，中是多粒子體系的勢函數(shù)，是中與體積V相應(yīng)的共扼動(dòng)量，是基本單元的體積，且由給出，為相應(yīng)的質(zhì)量，為外部壓強(qiáng)，代表各向異性應(yīng)力張量。2. 溫度恒定體系接下來討論滿足溫度恒定的基本方法。恒溫體系比定壓體系要抽象，簡單地給出視覺化圖像是困難的。勢能在原子坐標(biāo)和動(dòng)量之

23、外引入了附加自由度，并由下式用把現(xiàn)實(shí)系統(tǒng)與具有能量確定的(微正則系統(tǒng))假想系統(tǒng)聯(lián)系起來。(6.14)式中，為現(xiàn)實(shí)系統(tǒng)的時(shí)間，為假想系統(tǒng)的時(shí)間。該假想(力學(xué))系統(tǒng)的哈密頓量可以寫成(6.15)式中，是在假想系統(tǒng)中原子的動(dòng)量，是對應(yīng)的共扼動(dòng)量，是由公式給出的假想系統(tǒng)的自由度數(shù)目，二是熱浴的溫度，為玻爾茲曼常數(shù)。顯然，式(6.18)中右邊第一項(xiàng)代表粒子的動(dòng)能，第二項(xiàng)為粒子系統(tǒng)的勢能，第三項(xiàng)是熱浴的假想動(dòng)能，第四項(xiàng)是熱浴的勢能；因?yàn)樵诩傧胂到y(tǒng)中，力學(xué)系綜是微正則分布的，所以其配分函數(shù)變成具有總能量 E 的狀態(tài)總數(shù)(6.16)(6.17)若將此配分函數(shù)積分變換到現(xiàn)實(shí)系統(tǒng)中，則可得到在恒定溫度T下的正則分

24、布(也稱為玻爾茲曼分布)的配分函數(shù) (6.18)式中，是由下式給出的現(xiàn)實(shí)系統(tǒng)的哈密頓量(6.19)是假想系統(tǒng)的守恒量，而對真實(shí)系統(tǒng)來講，由于存在熱交換將不是守恒量。為了判斷計(jì)算結(jié)果是否合理，只要簡單地檢驗(yàn)是否保持恒定(守恒)就行了。3. NTP系綜到此，已經(jīng)導(dǎo)出了壓力恒定系綜和溫度恒定系綜的哈密頓量，將兩個(gè)系綜的哈密頓量組合起來就可以得到NTP系綜的哈密頓量，其模型如圖6-7所示。NTP系綜的哈密頓量為 (6.20)圖6-7 NTP系統(tǒng)示意圖若根據(jù)此哈密頓量導(dǎo)出正則方程式，并用速度置換動(dòng)量，則可得到各變數(shù)的運(yùn)動(dòng)方程(6.21)(6.22)(6.23)式中，變量上方的圓點(diǎn)表示對時(shí)間的微分，而分別

25、由下列各式給出(6.24)(6.25)(6.26)4. NPT系綜分子動(dòng)力學(xué)求解若在三維周期邊界條件下求解上述方程，則可模擬固體的結(jié)構(gòu)相變、溶解現(xiàn)象和晶化過程。各變量的運(yùn)動(dòng)方程式(6.21)，式(6.22)和式(6.23)是非線性常微分方程，不能用貝魯勒(Berrele)法進(jìn)行求解。因此，對這些方程通常用牛頓迭代方法求解。而問題的關(guān)鍵在于確定式(6.22) 和式(6.23)中的假想質(zhì)量。所考察基本單元內(nèi)的粒子體系存在著由粒子作用勢以及粒子數(shù)與粒子質(zhì)量共同決定的固有壓力和溫度等的起伏變化。若想再現(xiàn)這些波動(dòng)周期那樣給出和，則由體積變化或溫度變化引起的粒子系統(tǒng)對平衡態(tài)的偏離將迅速回復(fù)而變成穩(wěn)定狀態(tài)。

26、因此，若設(shè)此周期分別為和，則和可以表示為(6.27)(6.28)式中，為溫度，為玻爾茲曼常數(shù)，為原子數(shù)，為基本單元的線尺度，為壓縮率，和分別根據(jù)對NEV系綜的計(jì)算結(jié)果而定?！揪毩?xí)與思考】6-2. 推演NEV系綜的牛頓方程，哈密頓方程，說明各個(gè)物理量的意義，在那些實(shí)際的計(jì)算體系中得到應(yīng)用。6-3. 根據(jù)NPT系綜的推導(dǎo)方法，對NPH系綜的哈密頓方程做簡單的推導(dǎo)。6-4. 詳細(xì)的分析不同的系綜的使用范疇，給出各個(gè)系綜的使用條件，列表進(jìn)行簡單的比較。6.4 原子勢函數(shù)和分子力場構(gòu)造分子動(dòng)力學(xué)中存在著兩個(gè)基本的概念，一是在前面提到的系綜的問題，它所確定的是模擬系統(tǒng)的大小，規(guī)格和相關(guān)的模擬條件等信息；二

27、是本節(jié)需要講述的勢函數(shù)。簡單的講，勢函數(shù)就是用來描述原子和原子之間相互作用的。隨著勢函數(shù)的不斷發(fā)展，現(xiàn)有的勢函數(shù)主要分成了兩個(gè)大的部分。即是經(jīng)典的勢函數(shù)和電子理論的勢函數(shù)。如圖6-8所示，經(jīng)典分子動(dòng)力學(xué)的基本的勢函數(shù)的一個(gè)簡單分類。圖6-8 經(jīng)典的分子動(dòng)力學(xué)相互作用勢分類作用勢的選擇與動(dòng)力學(xué)計(jì)算的關(guān)系極為密切，選擇不同的作用勢，體系的勢能面會(huì)有不同的形狀，動(dòng)力學(xué)計(jì)算所得的分子運(yùn)動(dòng)和分子內(nèi)部運(yùn)動(dòng)的軌跡也會(huì)不同，進(jìn)而影響到抽樣的結(jié)果和抽樣結(jié)果的勢能計(jì)算，在計(jì)算宏觀體積和微觀成分關(guān)系的時(shí)候主要采用剛球模型的二體勢，計(jì)算系統(tǒng)能量，熵等關(guān)系時(shí)早期多采用Lennard-Jones、Morse勢等二體勢模型

28、，對于金屬計(jì)算，主要采用Morse勢，但是由于通過實(shí)驗(yàn)擬合的對勢容易導(dǎo)致柯西關(guān)系，與實(shí)驗(yàn)不符，因此在后來的模擬中有人提出采用EAM等多體勢模型，或者采用第一性原理計(jì)算結(jié)果通過一定的物理方法來擬合二體勢函數(shù)。但是相對于二體勢模型，多體勢往往缺乏明確的表達(dá)式，參量很多，模擬收斂速度很慢，給應(yīng)用帶來很大的困難，因此在一般應(yīng)用中，通過第一性原理計(jì)算結(jié)果擬合勢函數(shù)的Lennard-Jones，Morse等勢模型的應(yīng)用仍然非常廣泛。6.4.1 經(jīng)典理論的原子勢函數(shù)在分子動(dòng)力學(xué)的框架之下，早期的勢函數(shù)是主要以描述原子和原子之間的相互作用而產(chǎn)生和發(fā)展起來的。早在1903年的時(shí)候，Mie G就研究了兩個(gè)粒子之間

29、的相互作用勢，他給出的勢函數(shù)是由兩項(xiàng)構(gòu)成的，一項(xiàng)代表原子之間的相互吸引，另一項(xiàng)則代表著原子之間的相互排斥。這就是非常著名的Lennard-Jones勢函數(shù)的原始思考。對于其他有關(guān)勢函數(shù)的發(fā)展在這里就不做過多的贅述了。從圖6-9可以清楚的知道原子之間的相互作用可以分為四類：對勢、對泛函勢、組合勢、組合泛函勢。下面就對相關(guān)的勢函數(shù)的基本的發(fā)展做一個(gè)簡單的總結(jié)。表6-1 原子中相互作用勢函數(shù)勢函數(shù)作用因素典型代表對勢僅僅決定于原子的坐標(biāo)L-J勢、BMH勢、GK勢對泛函勢多體相互作用EAM勢有效介質(zhì)理論組合勢共價(jià)鍵的相互作用SW勢組合泛函勢Abell-Tersoff勢CheliKowsky-Phill

30、ips勢6.4.2 分子力場分子力場根據(jù)量子力學(xué)的波恩奧本海默近似，一個(gè)分子的能量可以近似看作構(gòu)成分子的各個(gè)原子的空間坐標(biāo)的函數(shù)，簡單地講就是分子的能量隨分子構(gòu)型的變化而變化，而描述這種分子能量和分子結(jié)構(gòu)之間關(guān)系的就是分子力場函數(shù)。分子力場函數(shù)為來自實(shí)驗(yàn)結(jié)果的經(jīng)驗(yàn)公式，可以講對分子能量的模擬比較粗糙，但是相比于精確的量子力學(xué)從頭計(jì)算方法，分子力場方法的計(jì)算量要小數(shù)十倍，而且在適當(dāng)?shù)姆秶鷥?nèi)，分子力場方法的計(jì)算精度與量子化學(xué)計(jì)算相差無幾，因此對大分子復(fù)雜體系而言，分子力場方法是一套行之有效的方法。以分子力場為基礎(chǔ)的分子力學(xué)計(jì)算方法在分子動(dòng)力學(xué)、蒙特卡羅方法、分子對接等分子模擬方法中有著廣泛的應(yīng)用。

31、 1. 力場主要討論的對象及其構(gòu)成：以下參數(shù)構(gòu)成一套力場函數(shù)體系需要有一套聯(lián)系分子能量和構(gòu)型的函數(shù)，還需要給出各種不同原子在不同成鍵狀況下的物理參數(shù)，比如正常的鍵長、鍵角、二面角等，這些力場參數(shù)多來自實(shí)驗(yàn)或者量子化學(xué)計(jì)算。 v 鍵伸縮能：構(gòu)成分子的各個(gè)化學(xué)鍵在鍵軸方向上的伸縮運(yùn)動(dòng)所引起的能量變化；v 鍵角彎曲能：鍵角變化引起的分子能量變化；v 二面角扭曲能：單鍵旋轉(zhuǎn)引起分子骨架扭曲所產(chǎn)生的能量變化； v 非鍵相互作用：包括范德華力、靜電相互作用等與能量有關(guān)的非鍵相互作用； v 交叉能量項(xiàng)：上述作用之間耦合引起的能量變化。2. 常用力場函數(shù)和分類不同的分子力場會(huì)選取不同的函數(shù)形式來描述上述能量與

32、體系構(gòu)型之間的關(guān)系。到目前，不同的研究小組設(shè)計(jì)了很多適用于不同體系的力場函數(shù)，根據(jù)他們選擇的函數(shù)和力場參數(shù)，可以分為以下幾類：力場CHARMM力場CVFF力場MMX力場AMBER力場COMPASS力場MMF94力場CFF力場UFF力場Dreiding力場ESFF力場第一代力場第二代力場第三代力場圖6-9 分子力場的發(fā)展和分類v 傳統(tǒng)力場(第一代力場)a) AMBER力場：由Kollman課題組開發(fā)的力場，是目前使用比較廣泛的一種力場，適合處理生物大分子；b) CHARMM力場：由Karplus課題組開發(fā)，對小分子體系到溶劑化的大分子體系都有很好的擬合；c) CVFF力場：CVFF力場是一個(gè)可以

33、用于無機(jī)體系計(jì)算的力場；d) MMX力場：MMX力場包括MM2和MM3，是目前應(yīng)用最為廣泛的一種力場，主要針對有機(jī)小分子。v 第二代力場 第二代的勢能函數(shù)形式比傳統(tǒng)力場要更加復(fù)雜，涉及的力場參數(shù)更多，計(jì)算量也更大，當(dāng)然也相應(yīng)地更加準(zhǔn)確；a) CFF力場CFF力場是一個(gè)力場家族，包括了CFF91、PCFF、CFF95等很多力場，可以進(jìn)行從有機(jī)小分子、生物大分子到分子篩等諸多體系的計(jì)算；b) COMPASS力場由MSI公司開發(fā)的力場，擅長進(jìn)行高分子體系的計(jì)算；c) MMF94力場Hagler開發(fā)的力場，是目前最準(zhǔn)確的力場之一。v 通用力場(第一代力場)通用力場也叫基于規(guī)則的力場，它所應(yīng)用的力場參數(shù)

34、是基于原子性質(zhì)計(jì)算所得，用戶可以通過自主設(shè)定一系列分子作為訓(xùn)練集來生成合用的力場參數(shù)；a) ESFF力場MSI公司開發(fā)的力場，可以進(jìn)行有機(jī)、無機(jī)分子的計(jì)算；b) UFF力場可以計(jì)算周期表上所有元素的參數(shù)；c) Dreiding力場適用于有機(jī)小分子、大分子、主族元素的計(jì)算。3. 分子力場的計(jì)算對于分子之間的相互作用通常需要將計(jì)算體系的所有粒子納入考慮范疇之內(nèi)，那么分子之間的相互作用可以作為長程作用來處理。在長程相互作用中又分為周期性長程相互作用力和非周期性相互作用力。周期性的長程作用力通?？梢岳脽o網(wǎng)格算法，包含Ewald求和算法等；而對于非周期性的相互作用通常使用快速多極子方法計(jì)算。下面針對這

35、兩個(gè)最為常見的算法做簡單的敘述。v Ewald求和算法Ewald求和算法是Ewald于1921年提出，最初是以計(jì)算離子晶體的能量。此方法選定一個(gè)計(jì)算系統(tǒng)的盒子，盒子中的粒子除了和盒子內(nèi)部的粒子相互作用之外，還會(huì)與其鏡像系統(tǒng)中的粒子發(fā)生作用，并且鏡像系統(tǒng)和盒子內(nèi)部的系統(tǒng)完全相同。圖顯示的是計(jì)算的相互作用的系統(tǒng)，中央的黑色的區(qū)域表示的是計(jì)算的區(qū)域，周圍的鏡像隨著系統(tǒng)的距離逐漸變淡，那么這樣的一個(gè)極限作用區(qū)域就像一個(gè)球，這個(gè)球稱為Ewald球。假設(shè)作用的系統(tǒng)的邊長是，含有個(gè)相互作用的粒子，其鏡像的作用位置可以表示為()；其中都是正整數(shù)。v 快速多極子法(The Fast Multipole Meth

36、od)在一個(gè)開放的系統(tǒng)中，不是對所謂的晶格進(jìn)行求和，而是對要求計(jì)算的整個(gè)體系的粒子對進(jìn)行求和。那么這樣的一個(gè)系統(tǒng)不同位置粒子的靜電能可以簡單的表述為(6.29)因此這個(gè)體系的靜電作用和前文所講述的Lennard-Jones作用項(xiàng)具有相同的形式?？焖俣鄻O子法將需要計(jì)算的體系分為多個(gè)體積相同的格胞，由格胞內(nèi)部的所有的原子計(jì)算每一個(gè)格胞的多極(電荷，偶極，四極)。而格胞和格胞之間的相互作用是利用中央多極展開的方法進(jìn)行計(jì)算。因此中央多極展開的方法的條件是格胞之間的距離必須是大于一定距離的，因?yàn)榇笥谝粋€(gè)格胞以上的以多極展開的方式計(jì)算相互作用的力，而格胞內(nèi)部的原子對是以成對原子的作用力來計(jì)算的?！揪毩?xí)與思

37、考】6-5. 查找文獻(xiàn)，推導(dǎo)Lennard-Jones勢函數(shù)的公式。6-6. 推導(dǎo)EAM勢函數(shù)的基本公式，并且實(shí)現(xiàn)簡單的算法，編寫基本設(shè)計(jì)程序。6-7. 以對勢中的Lennard-Jones勢函數(shù)為例，寫出對勢中勢函數(shù)的基本算法，并且進(jìn)行基本的程序設(shè)計(jì)。6.5 數(shù)值方法與編程技巧6.5.1 邊界條件隨著材料的低維化和納米化，材料的表面和邊界對于其性質(zhì)的影響愈來愈大。在處理物質(zhì)表面、界面，以及塊體物質(zhì)時(shí)，對其周期邊界條件問題的考慮就變得非常必要了。對于周期性邊界條件可以分為一維、二維及三維的情況。介紹薄膜(可使用二維周期邊界條件)和擴(kuò)大到半無限的表面情況，然后討論液體與晶體、晶體與晶體之間形成的

38、界面等，最后學(xué)習(xí)塊狀物質(zhì)的三維周期邊界條件。1. 薄膜的邊界條件(二維周期邊界條件)對于薄膜情況，如圖6-10所示賦予其周期性邊界條件?？梢哉J(rèn)為薄膜在x-y平面內(nèi)無限擴(kuò)展，存在周期邊界條件，而在 z方向受到限制(不賦予周期邊界條件)。事實(shí)上，二維周期性邊界條件也適用于在z方向不具有真正意義上的自由度的二維物質(zhì)。但就物質(zhì)科學(xué)的研究范疇而言，多數(shù)情況下是指在某種襯底上生成的原子或分子的薄膜。此處討論的邊界條件，除針對薄膜以外，也可近似地用于下面將要講到的在 z 的負(fù)方向半無限地?cái)U(kuò)展的表面問題，例如表面重構(gòu)等。但是，當(dāng)在z方向存在異質(zhì)原子層時(shí)，原來塊體物質(zhì)在異質(zhì)原子層附近將有近于表面的性質(zhì)，這種由塊

39、體向表面的過渡特征取決于原子層的厚度，一般要求原子層厚度在力的中斷距離以上(作用力程以上)。在處理實(shí)際問題時(shí)，只要嘗試進(jìn)行幾次改變原子層厚度的模擬和計(jì)算，并檢驗(yàn)計(jì)算結(jié)果是否變化就行了。圖6-10 薄膜周期性邊界條件示意圖2. 表面邊界條件(半無限的情況)對于表面存在表面重構(gòu)的情況，是與薄膜的情況不同的，其在 z 軸方向擴(kuò)展至半無限。處理此類系統(tǒng)之邊界條件時(shí)，只對x-y平面和 z 的負(fù)方向賦予相應(yīng)的周期性邊界條件，其示意圖如圖6-11所示。但是，在將基本單元原樣地復(fù)制中，必須考慮到這樣的一種情況，那就是在距薄膜結(jié)構(gòu)的基本單元最下層的相鄰底層被配置了基本單元表面層的復(fù)制品，從而導(dǎo)致了在本來應(yīng)具有塊

40、體物質(zhì)性質(zhì)的基本單元的最下層的結(jié)構(gòu)中出現(xiàn)了表面層的強(qiáng)烈影響。所以將基本單元對z軸反轉(zhuǎn)(即鏡像映射)進(jìn)行復(fù)制后使用。然而，即使在這種情況下，考慮到基本單元在z方向的厚度對計(jì)算結(jié)果有影響，所以應(yīng)當(dāng)進(jìn)行與前面所述同樣的檢驗(yàn)。圖6-11 采用映射方法半無限表面的周期性邊界條件圖6-12 固定層法半無限表面的周期性邊界條件在利用分子動(dòng)力學(xué)方法對分子束外延的晶體生長模擬時(shí)，要充分認(rèn)識到分子束外延 (MBE)是一個(gè)原子數(shù)不斷增加的非平衡系。所謂MBE技術(shù)就由束源產(chǎn)生的原子蒸發(fā)到襯底上，從而形成晶體。因?yàn)橛墒闯鰜淼脑泳哂袆?dòng)量，所以碰撞襯底，襯底的重心就變得下移，為了克服這種現(xiàn)象的影響，對分子束外延即在z軸

41、方向不賦予周期性邊界條件，而采用通過調(diào)整配置使原子位置坐標(biāo)不變的所謂固定原子層方法，然后使用邊界條件，如圖6-12所示。固定(原子)層還可理解為具有無限大質(zhì)量的原子層。來自固定層中原子作用力的貢獻(xiàn)可作為外力處理。為了基本單元最下層具有塊體物質(zhì)的性質(zhì)，固定層中的原子層數(shù)目必須足夠多。采用鏡像映射的方法可自然地處理體系伴隨溫度變化而產(chǎn)生的膨脹和收縮。但就使用固定層的方法而言，基本單元發(fā)生變化的同時(shí)，必須將固定層的晶格常數(shù)作為外部參量變化來考慮?！揪毩?xí)與思考】6-8. 通過對比映射方法和固定層法的不同，理解在分子動(dòng)力學(xué)中的二維邊界條件的設(shè)置方法，思考在線形的一維約束條件的邊界的處理方法。6-9. 根

42、據(jù)界面的邊界條件的設(shè)置方法，分析討論在分子動(dòng)力學(xué)中如何構(gòu)造線缺陷模型，做要的圖示和語言表述，說明界面的條件的選取方法。6-10. 根據(jù)上述的邊界的處理方法，試思考塊體材料的邊界條件的實(shí)現(xiàn)方法，并且給出簡單的實(shí)例。6.5.2 時(shí)間積分處理計(jì)算機(jī)模擬方法的基點(diǎn)是利用現(xiàn)代計(jì)算機(jī)高速和精確的優(yōu)點(diǎn)，對幾百個(gè)以至上千分子的運(yùn)動(dòng)方程進(jìn)行數(shù)值積分。有許多不同的積分方法，它們的效率和方便程度各異。基本問題就是用有限差分法來對二階常微分方程進(jìn)行積分。而對于分子動(dòng)力學(xué)而言其基本的工作也就是對運(yùn)動(dòng)方程進(jìn)行時(shí)間上的積分。這樣需要將運(yùn)動(dòng)方程離散化為有限的差分方程，即是將時(shí)間離散化為有限的格點(diǎn)。主要的方法是利用泰勒展開采取

43、不同的截?cái)喟霃?。而相鄰的格點(diǎn)之間的距離就是所謂的時(shí)間步長，這樣在知道了原始的位置坐標(biāo)之后，在知道某些物理量對于時(shí)間t的微分量之后就可以很容易確定下一個(gè)時(shí)間體系的坐標(biāo)。這里很容易出現(xiàn)所謂的兩類誤差，即是截?cái)嗾`差和舍入誤差。為此我們需要在實(shí)際的模擬中，應(yīng)該根據(jù)不同的問題來確定所需要的精度，為此發(fā)展了幾種時(shí)間積分的方法。1. Verlet方法Verlet法經(jīng)常用于NEV系綜的分子動(dòng)力學(xué)微分方程的數(shù)值積分求解。若將在和的位置坐標(biāo)分別用時(shí)刻t的位置坐標(biāo)進(jìn)行泰勒展開，則有(6.30)(6.31)為簡單起見，以后均省略有關(guān)原子標(biāo)識符號I，V(t)為速度，為了消去速度項(xiàng)，將式(6.30)和式(6.31)相加，

44、并將r()移至右邊，同時(shí)利用，則有(6.32)上式(6.32)是差分公式，其中()的奇次項(xiàng)全被抵消，同時(shí)忽略及其他高次項(xiàng)。差分精度(誤差)為0()。若由式(6-36)減去式(6-37)則可得到速度表達(dá)式(6.33)因?yàn)槭窍チ岁P(guān)于的偶數(shù)方項(xiàng)，所以速度的誤差為的量級。上述式(6.32)右邊第一、二項(xiàng)均是的量級項(xiàng)，而第三項(xiàng)是量級項(xiàng)。這樣若將大的絕對值項(xiàng)和小的絕對項(xiàng)混合，則隨著重復(fù)計(jì)算次數(shù)的增多，其誤差發(fā)生積累。為防止誤差積累，Hockney提出了所謂蛙躍(Leap-frog)法。2. 蛙躍法首先，利用速度將二階常微分方程轉(zhuǎn)換為二個(gè)一階常微分方程。亦即(6.34)速度的微分由時(shí)刻，和的差分給出(6.

45、35)根據(jù)此式，在時(shí)刻的速度給出如下(6.36)另一方面，利用時(shí)刻和時(shí)刻的差分可給出位置矢量的微分(6.37)根據(jù)上式可給出時(shí)刻的位置矢量表達(dá)式(6.38)上式(6.38)所表示的意義是速度和位置矢量的變化依賴步長。利用時(shí)刻的位置坐標(biāo)可計(jì)算出力，利用此力的結(jié)果求出時(shí)刻()的速度，進(jìn)而，利用此速度求出在()時(shí)刻的位置坐標(biāo)，這樣就完成了第一步的計(jì)算。正如以上所述，如果使用Verlet法或蛙躍法，則可求解關(guān)于NEV系綜的線性常微分方程。但是，在此之外的NTV，NHP，NTP等系綜給出的是非線性常微分方程。所以用Verlet法或蛙躍法不能進(jìn)行求解。關(guān)于這些非線性常微分方程的求解要采用下面介紹的所謂Ge

46、ar方法。Hockey提出的Leap-frog算法是Verlet算法的變化，這種方法涉及半時(shí)間間隔的速度。這種算法與Verlet算法相比有兩個(gè)優(yōu)點(diǎn)：(1)包括顯速度項(xiàng)；(2)收斂速度快，計(jì)算量小。這種算法明顯的缺陷是位置和速度不同步。3. Gear方法因?yàn)樵榆壽E是時(shí)間的連續(xù)函數(shù)，所以可將()時(shí)刻的位置、速度、加速度等對時(shí)刻泰勒展開(6.39)若使用給出的新的位置矢量，則可計(jì)算()時(shí)刻的力F()和加速度()。根據(jù)加速度可評價(jià)預(yù)測的加速度的誤差，亦即(6.40)若將此差值加于預(yù)測項(xiàng)上，則可得到修正項(xiàng)(6.41)(6.42)4. Velocity-Verlet算法這種算法可以同時(shí)給出位置、速度和加

47、速度，并且不犧牲精度，給出了顯速度項(xiàng)，計(jì)算量適中，目前應(yīng)用比較廣泛。(6.43)5. Gear的預(yù)測-校正算法這種算法分三步來完成。首先，根據(jù)泰勒展開，預(yù)測新的位置、速度和加速度。然后，根據(jù)新的計(jì)算的力計(jì)算加速度。最后這個(gè)加速度再由與泰勒級數(shù)展開式中的加速度進(jìn)行比較，兩者之差在校正步里用來校正位置和速度項(xiàng)。這種方法的缺點(diǎn)就是占有計(jì)算機(jī)的內(nèi)存大?！揪毩?xí)與思考】6-11. 查找文獻(xiàn)，總結(jié)分子動(dòng)力學(xué)中現(xiàn)有的時(shí)間算法，列表說明各種不同的時(shí)間算法優(yōu)勢。6-12. 推導(dǎo)出Verlet算法的基本方程，并編程實(shí)現(xiàn)。6-13. Hockey提出的Leap-frog算法是Verlet算法的變化，這種方法涉及半時(shí)間

48、間隔的速度，查找原始文獻(xiàn)，推導(dǎo)出Leap-frog算法的基本方程，結(jié)合練習(xí)6-12中的問題，修改程序?qū)崿F(xiàn)Leap-frog算法。6-14. 根據(jù)文中的敘述，詳細(xì)的推導(dǎo)出Gear算法的基本方程，并編程實(shí)現(xiàn)。6-15. 查閱文獻(xiàn)，分析討論時(shí)間可逆算法劉維公式在分子動(dòng)力學(xué)中應(yīng)用。6.6 計(jì)算結(jié)果的解析方法總體來說分子動(dòng)力學(xué)(MD)目前能夠討論的問題已經(jīng)可以涉及到比較廣闊的范疇，從基本的模擬體系的位形，到現(xiàn)在的各種材料的光譜學(xué)性能的計(jì)算。但是總體來說MD的計(jì)算結(jié)果僅僅分了兩個(gè)層次。第一個(gè)是計(jì)算模擬直接的輸出，基本上包含了在平衡態(tài)的體系位形信息和體系運(yùn)動(dòng)的基本動(dòng)力學(xué)信息；第二個(gè)是基于現(xiàn)代統(tǒng)計(jì)熱力學(xué)的平均

49、統(tǒng)計(jì)得到的一些宏觀體系的綜合信息，包含了光學(xué)，電學(xué)，磁學(xué)和光譜學(xué)等等。但是，在材料科學(xué)中，不僅要很好地知道原子水平上的三維結(jié)構(gòu)，而且還必須計(jì)算出表征宏觀物性(熱力學(xué)性質(zhì)、光學(xué)性質(zhì)等)的物理量值。為了靈活應(yīng)用分子動(dòng)力學(xué)方法，下面對與這些物性相關(guān)的物理參數(shù)值的計(jì)算方法作介紹。6.6.1 宏觀參量的統(tǒng)計(jì)平均由分子動(dòng)力學(xué)基本方程出發(fā)進(jìn)行計(jì)算的輸出結(jié)果是關(guān)于原子的位置坐標(biāo)和速度。物性參數(shù)值可根據(jù)位置坐標(biāo)和速度通過統(tǒng)計(jì)處理求出。在統(tǒng)計(jì)物理中，可利用系綜微觀量的統(tǒng)計(jì)平均值來計(jì)算物性參數(shù)值，亦即： (6.44)式中，ens是系綜(Ensemble)的英文縮寫，是分布函數(shù)(幾率密度)。幾率密度因所考察的系綜不同

50、而不同。例如，對于微正則系綜(NEV系)，可表達(dá)為(6.45)而對于正則系綜(NTV系)，等溫、等壓系綜(NPT系)和巨正則系綜，分別由下列各式給出：(6.46)系綜之間的演變可由相應(yīng)的公式聯(lián)系起來。首先，配分函數(shù)Z和熱力學(xué)勢，分別由下列式子給出(6.47)式中，在(,)中，設(shè)定，則可有以下組合：NTV和NEV系綜之間的轉(zhuǎn)變NTP和NTV系綜之間的轉(zhuǎn)變VT和NTV系綜之間的轉(zhuǎn)變在上面最后一個(gè)演變例子中，關(guān)于N的積分應(yīng)理解為對N0，1，的求和。在NEV系綜情況下，表示所有狀態(tài)數(shù)W，熱力學(xué)勢為(-TS)；對于NTP和VT系綜，其熱力學(xué)勢分別為N和(-PV)。此時(shí)，物理量在系綜中的平均值的變換可由下

51、式給出(6.48)在熱力學(xué)極限情況下，上述各系綜的平均是相等的(6.49)到此，容易發(fā)現(xiàn)，有了上述這些分布函數(shù)，代入式(6.50)進(jìn)行積分運(yùn)算，理論上講，可以求出各種各樣的物理量。然而，在通常情況下，式(6.50)不能給出真正的解析解。在理想的剛體球勢情況下，可以進(jìn)行解析求解，而對于連續(xù)函數(shù)勢，只能是數(shù)值求解。不難看出，這些解始終都是近似解。從歷史的角度來說，正是為了超越上述積分方程帶來的限制，人們引人了所謂蒙特卡羅方法和分子動(dòng)力學(xué)方法。在分子動(dòng)力學(xué)方法中，人們使用了所謂的時(shí)間平均等于系綜平均的各態(tài)歷經(jīng)假設(shè)，亦即(6.49)在熱力學(xué)統(tǒng)計(jì)物理學(xué)中，各態(tài)歷經(jīng)假設(shè)沒有嚴(yán)格地證明，它的正確性是由它的推

52、論被實(shí)驗(yàn)結(jié)果所證實(shí)而被肯定。式(6.59)中由下式給出 (6.50)在氣體動(dòng)力論中，原子的速度分布為Maxwell分布。在統(tǒng)計(jì)力學(xué)中，瞬間的熱力學(xué)性質(zhì)可以就其平均值做Gaussian分布描述。所以，對任意一個(gè)狀態(tài)函數(shù)X(t)，其瞬間的值是作無規(guī)地起伏變化，而此狀態(tài)值X(t)則依據(jù)它的平均值定義了其分布函數(shù)。但是要使得性質(zhì)分布函數(shù)X(t)的取樣點(diǎn)為可幸賴的，則模擬的向空間長度要多長？在此，先來分辨三種形式的統(tǒng)計(jì)平均值：軌跡平均(trajectory average)、取樣平均(sample average)、取樣平均的期望值(the expected value of the sample av

53、erage)。軌跡平均：考慮一個(gè)獨(dú)立系統(tǒng)的平衡相空間軌跡，并且想象在軌跡上取個(gè)不同的取樣點(diǎn)，這些取樣點(diǎn)是均勻分布在軌跡在線，相隔時(shí)間的相鄰點(diǎn)。由于軌跡是代表所有的平衡分布，而且很小，所以是一個(gè)非常大的值。對一個(gè)與時(shí)間有關(guān)的性質(zhì)物理量，時(shí)間平均值可取為軌跡平均值是我們想要知道的量。對于在任何時(shí)刻的相對于平均值的變異數(shù)為(6.51)取樣平均：雖然我們想要知道系統(tǒng)性質(zhì)的平均值(也就是軌跡平均值)，但是卻又不希望對所有軌跡上的做計(jì)算。所以找了一個(gè)取代的方法：在軌跡上取一小段，并在此段上取M點(diǎn)做計(jì)算(M)。也就是在一小段軌跡上取連續(xù)相鄰的M個(gè)點(diǎn)，使成一個(gè)集合，此集合稱為樣本。這個(gè)樣本有它自己的平均值(6

54、.52)而在軌跡上任意時(shí)刻的，可針對此平均值做變異數(shù)：(6.53)但是由軌跡中某段所得的平均值是否可代表整個(gè)軌跡的平均值？這個(gè)問題等介紹完第三個(gè)平均值后就可討論。取樣平均的期望值：自整個(gè)軌跡上取個(gè)片段軌跡，每段軌跡都包含M個(gè)點(diǎn)，任意其中一個(gè)取樣段(例如第k個(gè))，都有它自己的平均值，所以共有個(gè)平均值。取這些個(gè)平均值再做平均稱為平均值的期望值，則這些取樣段平均值會(huì)對做性質(zhì)函數(shù)分布，稱為取樣分布(sampling distribution)。對任意一個(gè)取樣段的平均值相對平均值的期望值的變異數(shù)為(6.54)由以上三種取樣(統(tǒng)計(jì))方法，可得到取樣理論的一個(gè)重要結(jié)果：期望值的變異數(shù)與軌跡變異數(shù)S成比例關(guān)系

55、(6.55)所以，若對整個(gè)軌跡取樣，則M=，因此，任一小段的平均值(近似期望值)是等于軌跡平均值。另外，若M，則(6.56)取任一段的平均值去估計(jì)軌跡平均值得準(zhǔn)確度是與任一小段上的取樣點(diǎn)數(shù)M有關(guān)，而與全軌跡長度的無關(guān)。直覺上，我們以為在取軌跡中的一段做平均值時(shí)，只有此段長度愈長(也就是)，其平均值才會(huì)愈接近全軌跡的平均值。但由上式卻看到，當(dāng)非常大時(shí)，是不重要的。這也就是為什么取1500人做調(diào)查的結(jié)果，可代表數(shù)千萬人的分布。將取樣理論應(yīng)用在分子動(dòng)力學(xué)：在分子動(dòng)力學(xué)中，常常針對時(shí)間及粒子數(shù)目做性質(zhì)函數(shù)的平均值。通常所謂的巨觀尺度，就時(shí)間而言，數(shù)量級為秒；就粒子數(shù)目而言，數(shù)量級為。但是在做分子動(dòng)力模

56、擬時(shí)，我們卻可以僅做皮秒及數(shù)百個(gè)分子的計(jì)算即可。【練習(xí)與思考】6-16. 搜索利用分子動(dòng)力學(xué)方法計(jì)算之后的位形顯示軟件，并且列表比較其優(yōu)勢。6-17. 根據(jù)文中提到的配分函數(shù)概念，利用熱力統(tǒng)計(jì)學(xué)的基礎(chǔ)知識，完成定壓比熱容CP，定容比熱容Cv，熱膨脹系數(shù)，等溫壓縮率參數(shù)的公式推到，并利用程序來實(shí)現(xiàn)。6-18. 根據(jù)文獻(xiàn)查找結(jié)果，編程完成一般分子動(dòng)力學(xué)方法中輸運(yùn)系數(shù)的計(jì)算。6-19. 根據(jù)擴(kuò)散系數(shù)的基本方程，通過查找參數(shù)模擬在Si體系中三個(gè)的H原子的擴(kuò)散的路徑，并編程實(shí)現(xiàn)。6.7 分子動(dòng)力學(xué)計(jì)算實(shí)例6.7.1 CNT的幾何特性自1991年被發(fā)現(xiàn)以來，碳納米管因其所特有的優(yōu)異物理、化學(xué)和機(jī)械性能及其巨大的潛在應(yīng)用價(jià)值而得到了廣泛的關(guān)注。例如根據(jù)碳納米管獨(dú)特的一維管狀結(jié)構(gòu)和極大的長徑比可將其制成高效的傳質(zhì)單元；超薄的針形尖端使碳納米管可用做掃描隧道顯微鏡及原子力顯微鏡的探針；碳納米管具有較強(qiáng)的場發(fā)射性能，是一種優(yōu)良的場效應(yīng)晶體管材料，Duan等最近在研究中發(fā)現(xiàn)，用金屬銫對單壁碳納米管進(jìn)行修飾能加強(qiáng)場發(fā)射性能；碳納米管兼具金屬或半導(dǎo)體導(dǎo)電性，是一種理想的電極材料；碳納米管的表面積較大，可應(yīng)用于儲(chǔ)氫、儲(chǔ)能以及吸附劑等領(lǐng)域。如上所述，碳納米管的發(fā)現(xiàn)進(jìn)一步拓寬了人們對納米材料的認(rèn)知領(lǐng)域，隨著制備與