分子模擬的一些基本資料

1、第1章 前 言1.1 研究背景1.1.1 分子模擬及其發(fā)展分子模擬(Molecular Simulation)為二十世紀(jì)初發(fā)展起來的一種計算機(jī)模擬方法，它泛指用于模擬分子或分子體系性質(zhì)的方法，主要用于探索研究具有三維結(jié)構(gòu)的分子結(jié)構(gòu)和分子的性能第2章 分子動力學(xué)模擬法分子模擬作為一種計算機(jī)模擬技術(shù)，主要可以進(jìn)行解釋工作和預(yù)測工作。前者為實驗奠定理論基礎(chǔ)，通過模擬解釋實驗現(xiàn)象、建立理論模型、探討過程機(jī)理等，后者為實驗過程提供可能性和可行性研究，進(jìn)行方案輔助設(shè)計、材料性能預(yù)測、過程優(yōu)化篩選等。分子模擬所涉及的研究領(lǐng)域，涵蓋了物理、化學(xué)、化工、材料、生化等幾乎一切可以通過建立理論模型進(jìn)行研究的體系，多

2、數(shù)能夠得到與實驗結(jié)果近似的計算結(jié)果，所以分子模擬己經(jīng)逐漸成為與實驗技術(shù)并重的強(qiáng)有力的研究手段。20世界90年代以來，由于計算機(jī)科學(xué)技術(shù)的飛速發(fā)展，分子模擬的地位日顯突出。現(xiàn)在，理論分析、實驗測定和模擬分析已成為現(xiàn)代科學(xué)的三種主要方法，在化工新產(chǎn)品、新材料和新藥物的研究和開發(fā)中，采用分子模擬技術(shù)，可以從分子的微觀性質(zhì)推算及預(yù)測產(chǎn)品及材料的介觀、宏觀性質(zhì)，已成為新興學(xué)術(shù)方向。為探求納米尺度下傳熱現(xiàn)象的規(guī)律和內(nèi)在機(jī)理，我們必須從微觀細(xì)節(jié)著手，研究載熱粒子的行為，并依據(jù)統(tǒng)計力學(xué)原理得到系統(tǒng)的統(tǒng)計性質(zhì)。分子模擬技術(shù)為某些問題(如統(tǒng)計力學(xué))的精確解答起到非常重要的作用，否則這些問題很難或只能近似求解。此外

3、，一些在極溫極壓條件下很難或者根本無法進(jìn)行的實驗，利用分子模擬卻比較容易做到，在這種情況下，它實際上是對理論應(yīng)用范圍的一個較好的檢測。通過對比分子模擬的結(jié)果與真實的實驗結(jié)果，分子模擬的結(jié)果可以檢驗理論的正確性，相反，真實的實驗結(jié)果與理論可以檢測模型構(gòu)建的合理性，此兩者相互聯(lián)系相互進(jìn)步，通過對比，我們可以修正所建立的模型，并從模擬結(jié)果中進(jìn)一步加深對所研究問題的理解，有助于解釋出現(xiàn)的新現(xiàn)象。2.1 分子動力學(xué)概述分子動力學(xué)模擬（Molecular Dynamics Simulation ,MDS）是目前最廣泛計算龐大復(fù)雜體系所采用的方法，并且現(xiàn)已廣泛應(yīng)用于計算物理、化學(xué)、材料、生物科學(xué)等領(lǐng)域內(nèi)的一

4、種計算機(jī)模擬技術(shù)，它的思想是將組成系統(tǒng)的微觀粒子視為牛頓經(jīng)典粒子，將所研究的系統(tǒng)視為經(jīng)典多體系統(tǒng)，通過選擇恰當(dāng)?shù)膱龊瘮?shù)來描述系統(tǒng)內(nèi)微觀粒子間相互作用的勢函數(shù)及系統(tǒng)外加約束條件后，利用牛頓運(yùn)動定律來求解微觀粒子的牛頓運(yùn)動方程。在MDS中，勢能模型的選取是非常重要的環(huán)節(jié)。模擬能否成功取決于能否正確地選擇勢能模型。對于簡單分子，常采用硬球、軟球、Lennard-Jones模型、Kihara模型等等。對復(fù)雜分子，可以采用多中心的位置位置相互作用模型，但中心間的相互作用仍為簡單模型，對帶電分子或離子，還要引入庫侖作用，常用的多體作用勢有Stillinger-Weber(SW)作用勢、Tersoff勢等。

5、自1966年起，由于分子力學(xué)的發(fā)展，人們又系統(tǒng)的建立了許多適用于生化分子體系、聚合物、金屬及非金屬材料的力場，使得計算復(fù)雜體系的結(jié)構(gòu)與一些熱力學(xué)、光譜性質(zhì)的能力及精準(zhǔn)性大為提升。分子動力學(xué)模擬是應(yīng)用這些力場及根據(jù)牛頓運(yùn)動力學(xué)原理所發(fā)展的計算方法。此方法的優(yōu)點為精準(zhǔn)性高，同時可以獲得系統(tǒng)的動態(tài)與熱力學(xué)統(tǒng)計資料，并可廣泛地適用于各種系統(tǒng)及各類特性的研究。與蒙地卡羅計算方法相比較，分子動力學(xué)模擬系統(tǒng)中粒子的運(yùn)動有正確的物理依據(jù)，并且計算技巧經(jīng)過許多修改，現(xiàn)已日趨成熟，由于其計算能力強(qiáng)，能滿足各類問題的需求，許多使用方便的分子動力學(xué)模擬商業(yè)軟件也已問世。2.2 分子動力學(xué)模擬的基本步驟分子動力學(xué)模擬（

6、Molecular Dynamics Simulation ,MDS）的出發(fā)點是將組成系統(tǒng)的微觀粒子視為牛頓經(jīng)典粒子，將所研究的系統(tǒng)視為經(jīng)典多體系統(tǒng)，通過選擇恰當(dāng)?shù)膱龊瘮?shù)來描述系統(tǒng)內(nèi)微觀粒子間相互作用的勢函數(shù)及系統(tǒng)外加約束條件后，利用牛頓運(yùn)動定律來求解微觀粒子的牛頓運(yùn)動方程，的到系統(tǒng)隨時間演進(jìn)的微觀過程，最后基于統(tǒng)計力學(xué)計算得到系統(tǒng)的各種參數(shù)。分子動力學(xué)模擬一般由下面幾個步驟組成：1.選取要研究的系統(tǒng)及其邊界和系統(tǒng)內(nèi)粒子間的作用勢（力場函數(shù)）；2.設(shè)定系統(tǒng)中粒子的初始位置和動量初始化；3.建立模擬算法，計算粒子間作用力及各粒子的速度和位置，并求解Newton運(yùn)動方程，得到下一時刻粒子的空間位置

7、和速度；4.當(dāng)體系達(dá)到平衡后，依據(jù)相關(guān)的統(tǒng)計公式，獲得各宏觀參數(shù)和輸運(yùn)性質(zhì)；2.2.1 選擇恰當(dāng)系統(tǒng)與邊界分子動力學(xué)模擬是在一定的系綜下進(jìn)行的，常見的系綜有微正則系綜、正則系綜、等溫等壓系綜和等溫等焓系綜等。在統(tǒng)計物理中，系綜（Ensemble）代表一群類似的體系的集合。對一類相同性質(zhì)的體系，其微觀狀態(tài)（比如每個粒子的位置和速度）仍然可以大不相同。實際上，對于一個宏觀體系，所有可能的微觀狀態(tài)數(shù)是天文數(shù)字。統(tǒng)計物理的一個基本假設(shè)（各態(tài)歷經(jīng)假設(shè)）是：對于一個處于平衡的體系，物理量的時間平均等于對應(yīng)系綜里所有體系進(jìn)行平均的結(jié)果。體系的平衡態(tài)的物理性質(zhì)可以對不同的微觀狀態(tài)求和來得到。系綜的概念是由約西

8、亞·威拉德·吉布斯（J. Willard Gibbs）在1878年提出的。在分子動力學(xué)模擬過程中，經(jīng)常采用的系綜有微正則系綜(Microcanonical Ensemble)、正則系綜(Canonical Ensemble)。正則系綜 (Canonical Ensemble)也稱為NVT系綜，它是統(tǒng)計力學(xué)中系綜的一種。它代表了許多具有相同溫度的體系的集合，該系統(tǒng)具有恒定的粒子數(shù)N、體積V和溫度T。正則系綜是最普遍應(yīng)用的系綜。通常，系綜內(nèi)每個體系的粒子數(shù)和體積都是相同的。但每個體系都可以和系綜內(nèi)其他體系交換能量。同時系綜里所有體系的能量總和，以及所有體系的總個數(shù)是固定的。但是

9、正則系綜的恒溫條件是不能自發(fā)實現(xiàn)的，因此在分子動力學(xué)模擬過程中經(jīng)常會采用一定的措施來保證系統(tǒng)溫度的恒定。在分子動力學(xué)模擬過程中，系統(tǒng)的溫度T又Boltzmann能量均分定理得出，即:T=i=1Nmivi23NkB （2. 1）式中，N代表區(qū)域中原子的數(shù)量，vi為粒子i的速度，kB=1.38×10-23 J·K-1。此公式成立的前提是系統(tǒng)區(qū)域內(nèi)達(dá)到熱力學(xué)平衡，若系統(tǒng)區(qū)域內(nèi)無法達(dá)到熱力學(xué)平衡，溫度的定義就失去了意義。對于NVT系綜，調(diào)節(jié)溫度恒定的有效方法是直接調(diào)節(jié)系統(tǒng)內(nèi)微觀粒子的速度。在這些條件下，當(dāng)系綜內(nèi)所有體系被分配到不同的微觀態(tài)上，我們發(fā)現(xiàn)：每個微觀態(tài)上的體系個數(shù)正比于e

10、xp(-E)。其中=1/kBT，kB是玻爾茲曼常數(shù)，T是絕對溫度。正則系綜的配分函數(shù)是：Z=exp(-E)（2. 2）配分函數(shù)的對數(shù)就是亥姆霍茲自由能（Helmholtz Free Energy，符號F），當(dāng)配分函數(shù)計算出以后，平均能量可以直接從ln Z對一階導(dǎo)數(shù)中求得。F=-kBT lnZ（2. 3）<E>=-Z（2. 4）微正則系綜(Microcanonical Ensemble)又稱為NVE系綜，該系綜由許多具有相同能量E，粒子數(shù)N，體積V的體系的集合。微正則系綜可以通過平衡態(tài)分子動力學(xué)由系統(tǒng)內(nèi)微觀粒子的牛頓運(yùn)動方程自動實驗，比較簡單方便，因此多數(shù)平衡態(tài)分子動力學(xué)模擬是在微正

11、則系綜下進(jìn)行的。對于分子動力學(xué)模擬除了要選定合適的系綜之外，還需要對模擬體系加以適當(dāng)?shù)倪吔鐥l件。這是因為材料的宏觀性質(zhì)是由大量粒子的相互作用所表現(xiàn)出來的，在模擬過程中，如果所建立的模型含有足夠多的粒子，那么才能真實的模擬出宏觀材料的性質(zhì)。但是計算機(jī)的計算能力是有限的，分子動力學(xué)模擬時粒子數(shù)目一般不超過10000個，加上系綜的某些限制性條件（例如NVT系綜中，調(diào)節(jié)系綜溫度恒定的某些方法），模擬的粒子數(shù)目可能要比這個數(shù)還要小。這樣少的粒子所組成的系統(tǒng)與真實的宏觀材料體系是具有一定誤差的，同時粒子數(shù)目少，其表面原子與體內(nèi)原子之比顯然很大，因此會造成明顯的表明效應(yīng)。因此，為了減小小體系的影響且使計算中

12、系統(tǒng)的密度維持恒定，在分子動力學(xué)模擬過程中通常會采用有效的邊界條件。常用于分子動力學(xué)模擬中的邊界條件有三種：周期性邊界條件、對稱性邊界條件和固壁邊界條件。在周期性邊界條件下，位于圖中央的盒子表示所模擬的系統(tǒng)，其周圍的盒子與模擬系統(tǒng)具有相同的排列及運(yùn)動，我們稱之為周期性鏡像（Periodic mirror image）系統(tǒng)，當(dāng)計算系統(tǒng)中任意微粒移出盒子后，則必有一粒子由相對方向移入，這樣的限制條件使得系統(tǒng)中的粒子數(shù)目維持恒定，密度不變，符合實際要求。在周期性邊界條件下，計算系統(tǒng)中分子間的作用力變得復(fù)雜，我們通常采取最近鏡像法（Nearest mirror image）來簡化問題。如圖 2- 1

13、周期性邊界示意圖。圖 2- 1 周期性邊界示意圖2.2.2 分子間作用力的描述力場分子力場（Forcefields）根據(jù)量子力學(xué)的波恩-奧本海默近似可以將力場看作是粒子勢能面的經(jīng)驗表達(dá)式，即將一個分子的勢能近似看作構(gòu)成分子的各個原子的空間坐標(biāo)的函數(shù)，簡單地講就是分子的勢能隨分子構(gòu)型的變化而變化，而描述這種分子勢能和分子結(jié)構(gòu)之間關(guān)系的就是分子力場函數(shù)。分子動力學(xué)模擬過程中計算系統(tǒng)內(nèi)粒子間的作用力主要依靠的就是力場給各個粒子間分配的作用力模型，分子動力學(xué)模擬的主要依據(jù)是經(jīng)典力學(xué)，那么分子動力學(xué)模擬結(jié)果的準(zhǔn)確性也主要取決于力場函數(shù)對粒子間相互作用力描述的正確與否。由于分子動力學(xué)模擬研究的廣泛性，分子

14、間相互作用模型并非完全相同，針對不同的分子類型、分子體系，分別有不同的力場函數(shù)可以使用。分子力學(xué)與分子動力學(xué)發(fā)展至今已有40多年，其計算的系統(tǒng)由原來最初的單原子分子系統(tǒng)延伸至多原子分子、聚合物分子、生化分子系統(tǒng)。計算時所需的力場也隨著系統(tǒng)的復(fù)雜度的增加而增加，在執(zhí)行分子動力學(xué)計算時，選擇合適的力場十分重要。常見的力場有CVFF、CFF91、CFF95、PCFF、AMBER、CHARMM、COMPASS和MM形態(tài)力場等。MM形態(tài)力場為Allinger等所發(fā)展，此力場適合用于模擬有機(jī)化合物、自由基、離子等隨著分子模擬的進(jìn)步MM形態(tài)力場也發(fā)展出MM2、MM3和MM4形態(tài)力場等，MM力場將一些常見的原

15、子細(xì)分（例如，將碳原子分為sp3、sp2、sp、酮基碳、碳自由基和碳陽離子等并且這些不同形態(tài)的碳原子具有不同的力場參數(shù)），并為這些原子補(bǔ)充更多的力場參數(shù)，以更精確詳細(xì)的來描述相同類型、不同電子結(jié)構(gòu)原子的受力。CVFF力場其全名為一致性價力場（Consistent valence force field）,為Dauber Osguthope等所發(fā)展。此力場以計算系統(tǒng)的結(jié)構(gòu)與結(jié)合能最為準(zhǔn)確，亦可以提供合理的構(gòu)型能和振動頻率，常用于計算氨基酸、水及含各種官能團(tuán)的分子體系，后經(jīng)過不斷改進(jìn)，CVFF力場已適用于計算各種多肽、蛋白質(zhì)與大量的有機(jī)分子。AMBER力場（Assisted model build

16、ing with energy minimization）為美國加州大學(xué)的Peter Kollman等所發(fā)展，此力場通常可以得到分子的幾何結(jié)構(gòu)、構(gòu)型能、振動頻率與溶劑化自由能（Solvation free energy）等數(shù)據(jù),所以常用于較小的蛋白質(zhì)、核酸、多糖等生化分子的模擬計算。CHARMM力場（Chemistry at Harvard macromolecular mechanics）為美國哈佛大學(xué)所發(fā)展，此力場可研究許多分子系統(tǒng)，例如小的有機(jī)分子、溶液、聚合物、生化分子等。此力場參數(shù)除來自計算結(jié)果與實驗值的比對外，還引用了大量的量子計算結(jié)果。應(yīng)用此力場能夠得到與實驗值相近的結(jié)構(gòu)、作用能

17、、構(gòu)型能、轉(zhuǎn)動能、振動頻率、自由能與許多和時間有關(guān)的物理量。隨著力場應(yīng)用研究的深入，隨后出現(xiàn)了第二代力場。第二代力場的形式遠(yuǎn)比經(jīng)典里場復(fù)雜，并需要大量的力場參數(shù)。其設(shè)計的目的是為了更加精確的計算分子的各種性質(zhì)、結(jié)構(gòu)、光譜、熱力學(xué)特性、晶體特性等，其力場參數(shù)除了引用大量的實驗數(shù)據(jù)之外，并參照了精確的量子計算結(jié)果，尤其適用于有機(jī)分子系統(tǒng)。CFF91、CFF95、COMPASS與PCFF都屬于第二代力場。一般情況下，分子的總能量為動能與勢能的和，而分子的勢能可表示為簡單的幾何坐標(biāo)的函數(shù)。復(fù)雜分子的總勢能一般可分為各類型勢能的總和，即由成鍵項（Ubonded）、交叉項(Ucross)和非鍵項(Unob

18、onded)三部分構(gòu)成，即：總勢能=成鍵項勢能+交叉項勢能+非鍵項勢能習(xí)慣上常用的數(shù)學(xué)表達(dá)式為：U=Ubonded+Ucross+Unobonded成鍵項（Bonded term）包括鍵伸縮項(Ub)、鍵角彎曲項(U)、二面角扭曲項(U)和離平面振動項(U)。交叉能量項(Cross energy term)主要指成鍵作用之間耦合引起的能量變化，包括鍵伸縮-鍵伸縮(Ubond-bond)、鍵伸縮-鍵角彎曲-鍵伸縮(Ubond-angle-bond)、鍵角彎曲-鍵角彎曲(Uangel-angle)、二面角扭轉(zhuǎn)項-鍵伸縮(Utorsion-bond)、二面角扭轉(zhuǎn)項-鍵角彎曲-鍵角彎曲(Utorsio

19、n-angle-angle)、鍵角彎曲-二面角扭轉(zhuǎn)項-鍵角彎曲(Uangle-torsion-angle)、離平面振動-二面角扭轉(zhuǎn)項-離平面振動(Uoop-torsion-oop)，一般情況下，交叉能量項對整體勢能影響較小,所以在一般的力場中并不含有這些項,只有在計算結(jié)果精度要求較高時才有交叉能量項。交叉項通常形式比較復(fù)雜,計算起來非常耗時,尤其是在計算原子數(shù)目較多的體系時。非鍵項（Nobonded term）包括范德瓦耳斯作用（Uvdw）、庫倫靜電作用（Uel）和氫鍵作用（UH）等。以下著重強(qiáng)調(diào)幾個重要的勢能項： (1) 非鍵結(jié)范德瓦耳斯勢能(Uvdw) 若A、B兩原子屬于同一分子但其間隔多

20、于兩個連接的化學(xué)鍵（如 ACCB），或者兩原子屬于不同的分子，雖無直接連接成鍵但因為距離較近，則此兩原子間有非鍵結(jié)范德瓦耳斯力（Nobonded van der waals force）。計算非鍵結(jié)作用時，通常將某個原子視為位于其原子核坐標(biāo)的一點，力場中最常見的非鍵結(jié)勢能形式為Lennard-Jones(L-J)勢能，這種勢能又稱為12-6勢能，其數(shù)學(xué)表達(dá)式為：Ur=4r12-(r)6（2. 5）式中與為勢能參數(shù)（因原子對而異），r為原子對間的距離，r-12項為排斥項，r-6項為吸引項，當(dāng)r很大時勢能趨近于零，這表明原子對距離很遠(yuǎn)時已無非鍵結(jié)作用，如圖 2- 2 L-J勢能曲線，表 2. 1部

21、分常見分子的LJ勢能參數(shù)。 圖 2- 2 L-J勢能曲線表 2. 1部分常見分子的LJ勢能參數(shù)原子或分子（/kb）/ Kr/ He10.412.602Ar141.63.350O2126.33.382CH4161.33.721C3H8268.54.992注: 1 =0.1nm 。 (2)鍵伸縮勢能(Ub) 在分子中相互鍵結(jié)的原子形成化學(xué)鍵，如CH，OH鍵等，化學(xué)鍵的鍵長并非恒定不變，而是在其平衡值附近呈小幅度的振動，這種振動所產(chǎn)生的勢能稱為鍵伸縮勢能或鍵伸縮勢能項（Bond stretching term）。鍵伸縮勢能項的一般形式為(將此形式視為簡諧振動)：Ub=12ikb(ri-ri0)2（2

22、. 6）式中kb為鍵伸縮的彈力常數(shù)，彈力常數(shù)越大，振動頻率越高，振動越快；ri和ri0分別表示第i個鍵的鍵長與第i個鍵的平衡鍵長。表 2. 2 MM2力場一些鍵的伸縮勢能參數(shù)鍵結(jié)平衡鍵長r/ 彈力常數(shù)kb/kcal/(mol·2)Csp3Csp31.523317Csp3Csp21.497317Csp2=Csp21.337690Csp3H1.113662Csp2=O1.438777Csp3Nsp31.438367注：1kcal=4.187 kJ, 1 =0.1nm 。有些力場（如MM3力場）為了增加計算的準(zhǔn)確性，除了簡諧振動項外增加了非簡諧振動項(除了二次的簡諧振動項外，還增加了三次與

23、四次的非間歇振動項)以更準(zhǔn)確的描述鍵伸縮勢能，如：Ub=i12kb(ri-ri0)2+k3(ri-ri0)3+k4(ri-ri0)4（2. 7） (3) 鍵角彎曲勢能(U) 若分子中連續(xù)鍵結(jié)的三原子形成一定的鍵角。與鍵的伸縮一樣，這些鍵角并非維持不變，而是在其平衡值附近呈小幅度的振蕩（一般視為簡諧振蕩），通過鍵角彎曲項（Angle bending term）可以描述這項勢能項，鍵角彎曲項的一般形式為：U=12ik(i-i0)2（2. 8）式中i及i0分別表示第i個鍵的鍵角及其平衡鍵角，k為鍵角彎曲的彈力常數(shù)。為了增加計算的精確性，我們可以像鍵伸縮勢能項一樣增加高次非簡諧項以提高計算的精確度。(

24、4) 二面角扭曲勢能(U) 若分子中連續(xù)鍵結(jié)四個原子并形成二面角（Dihedral angle），如丙烷CH3CH2CH3中的HCCC、HCCH等，一般分子中的二面角易于扭動，所以需要一個勢能項來描述這種勢能，我們稱之為二面角扭曲項（Torsion angle term），二面角扭曲項的一般形式為：U=12iV11+cos+V21-cos2+V31+cos3（2. 9）式中，為二面角角度，V1、V2和V3均為二面角扭曲項的彈力常數(shù)。(5) 離平面振動勢能(U) 在分子中，有些原子有共平面的傾向，如丙酮CH3C(O)CH3中，碳原子與氧原子的平衡位置位于共同平面，其他如苯環(huán)上的碳?xì)湓?、烯類中?/p>

25、碳?xì)湓拥染泄财矫娴膬A向，通常共平面的四個原子的中心原子會離開平面小幅度的振動，所以我們使用離平面振動項（Out-of-plane bending term）來描述這種勢能，離平面振動項的一般形式為：U=12ik 2(2. 10)式中，為離開平面振動的角度，k為離開平面振動項的彈力常數(shù)。(6) 庫倫靜電勢能（Uel） 離子或帶有部分電荷的分子，則這些帶電荷的粒子間因相對距離存在著靜電吸引或排斥作用，我們通過庫倫作用項（Columbic interaction term）來描述這種靜電作用，庫倫作用項的一般形式為：Uel=i,jqiqjDrij(2. 11)式中，qi和qj為第i個粒子與第j個

26、粒子所帶的電荷，rij為所帶電荷粒子間的距離，D為有效介電常數(shù)（effective dielectric constant）。2.2.3 Lennard-Jones(L-J)作用勢與Tersoff作用勢對于單壁碳納米的分子動力學(xué)模擬，由于單壁碳納米管模型中既有兩體相互作用也有三體相互作用，所以針對此模擬體系，采用了兩種作用勢，分別是Lennard-Jones作用勢與Tersoff作用勢。Lennard-Jones(L-J)作用勢（又稱L-J勢能函數(shù)或6-12勢能函數(shù)）是計算化學(xué)中用來模擬兩分子間作用勢能的一個函數(shù)，是一種兩體作用勢。最早由曼徹斯特大學(xué)的數(shù)學(xué)家John Lennard-Jones

27、于1931年提出。由于其解析形式簡單所以被廣泛使用，L-J作用勢的表達(dá)式為： Vrij=4rij12-(rij)6（2. 12）L-J作用勢相應(yīng)的兩體作用力為：Fr=-Vr=-ddrVrijr=4(1212rij13-66rij7)r（2. 13）式中為L-J作用勢的勢能阱的參數(shù)，是L-J作用勢的互相作用勢能正好為零時的兩體距離（也稱為平衡常數(shù)）,rij為粒子i、j間的距離。Tersoff作用勢是一種典型的三體作用勢，它已被廣泛用于描述碳族化合物的原子間的作用力。Tersoff作用勢的公式如下：ij=ijferijVR(rij)-bij*VA(rij)（2. 14）其中：VRrij=S-1ex

28、p-2S(rij-)（2. 15）VArij= SS-1exp-2/S(rij-)（2. 16）ferij=1,rij<R-D 12-12sin2(rij-R)D,R-D<rij<R+D0,rij>R+D （2. 17）bij*=bij+bji2,bij=1+an(k(i,j)ferik g(ijk)n-（2. 18）g=1+c2d2-c2d2+(h-cos)2（2. 19）式中，、為L-J勢能參數(shù)，ijk表示粒子i-j鍵和粒子i-k鍵之間的夾角，鍵序bij包含多體相互作用。公式中參數(shù)的選取如表：表 2. 3 Tersoff勢函數(shù)參數(shù)eVSnm-1ancdhRnmDnm

29、5.16441.4471.576919.6401.5724×10-70.727511/(2n)3.80494.3484-0.570580.1950.015通過Tersoff作用勢函數(shù)就可對勢函數(shù)求導(dǎo)來計算各個原子的受力情況了。2.2.4 公式的無量綱化在分子動力學(xué)模擬計算系統(tǒng)中原子或分子的運(yùn)動時，會涉及到大量的浮點運(yùn)算和指數(shù)運(yùn)算，如原子質(zhì)量以g為單位時，則原子質(zhì)量的量綱為10-22g；位置若以cm為單位時，其量綱為10-8cm；積分步程以秒（s）為單位時，其量綱為10-1310-16 s。這些量的量綱都非常小，在計算機(jī)計算時會增加大量的浮點運(yùn)算，可能在模擬中導(dǎo)致計算的誤差，因此在實際

30、執(zhí)行分子動力學(xué)模擬計算時都會對計算公式進(jìn)行無量鋼化處理（也可看作是簡化單位）以減少誤差。本文中涉及到的簡化單位轉(zhuǎn)換式有：能量：E*=E（2. 20）時間：t*=tm（2. 21）溫度：T*=kBT（2. 22）原子間距離：r*=r（2. 23）速率：v*=vm（2. 24）作用力：f*=f（2. 25）左上角用*標(biāo)記的量就是簡化單位后的參數(shù)，在分子動力學(xué)模擬計算中，各變量應(yīng)按照簡化單位后的數(shù)值進(jìn)行運(yùn)算。2.2.5 粒子的位置和動量初始化在分子動力學(xué)模擬過程開始前，必須先確定所有粒子的初始位置和速度，為了能使模擬系統(tǒng)盡快達(dá)到平衡并減少計算所需的時間，分子的初始速度應(yīng)接近真實情況。由于系統(tǒng)經(jīng)過一定

31、的模擬時間達(dá)平衡后，分子的速度的分布近似滿足Maxwell統(tǒng)計速度分布，所以采取Maxwell統(tǒng)計分布來賦予分子的初始速度是非常理想的，Maxwell統(tǒng)計分布如下：pvix=(mi2kBT)1/2exp-12mivix22kBT （2. 26）該式提供了一個質(zhì)量為mi的原子在溫度T時，在X方向的速度vix的高斯分布。隨著模擬時間的推進(jìn)，模擬的系統(tǒng)將遵循能量最小化的原則達(dá)到穩(wěn)定狀態(tài)。2.2.6 求解粒子的Newton運(yùn)動方程在分子動力學(xué)模擬中，必須求解式的牛頓運(yùn)動方程以計算速度與位置。而運(yùn)動方程的迭代方法及離散方法很多，其中以Verlet算法、Gear校正預(yù)測法算法（Predictor-Corr

32、ector method）最為常用。在Verlet方法中，將粒子的位置以泰勒公式展開，即：rt+t=rt+ddtrtt+12!d2dt2rt(t)2+(2. 27)將式中得t換成-t可得：rt-t=rt-ddtrtt+12!d2dt2rt(t)2+(2. 28)將式2.27與式2.28相加得：rt+t=-rt-t+2rt+d2dt2rt(t)2(2. 29)因為d2dt2rt=a(t),所以根據(jù)式2.29，可由t及t-t的位置預(yù)測t+t時的位置。將式1.7與式2.27相減的速度式：vt=drdt=12trt+t-r(t-t) (2. 30)此式表示時間t時的速度可由t+t與t-t的位置得到。V

33、erlet式的缺點在也于此，因為式中含有1/t的誤差項，為了矯正此缺點，Verlet發(fā)展了一種跳蛙方法(Leap frog method)的計算式,此方法計算速度與位置的數(shù)學(xué)表達(dá)式為：vit+12t=vit-12t+ait(2. 31)rit+t=rit+vit+12tt(2. 32)計算時假設(shè)已知vit-12t與rit，則由t時的位置rit計算質(zhì)點所受的力與加速度ai(t)。再依式（）來預(yù)測時間為t+12t時的速度vit+12t，并以此類推。時間t時的速度可由式（）得出：vit=12 vit+12t+vit-12t (2. 33)利用Verlet跳蛙方法計算僅需儲存rit與vit-12t兩種

34、數(shù)據(jù)，因此可以節(jié)省儲存空間，這種方法使用簡單且穩(wěn)定性、準(zhǔn)確性都較高。2.2.7 導(dǎo)熱系數(shù)的計算通過基于Green-Kubo線性響應(yīng)理論的平衡態(tài)分子動力學(xué)（EMD）方法，可知導(dǎo)熱系數(shù)的表達(dá)式為：k=13VkBT20+<J(t)J(0)>dt(2. 34)式中，kB波爾茲曼常數(shù)，T溫度，V模擬區(qū)域的體積，J為有效熱流。系統(tǒng)的總熱流表達(dá)式為：Jqt=ddtiritei(t)=ivitei(t)+iritdeidt(2. 35)其中，vi為原子i的速度，ei為原子所具有的能量。當(dāng)分子動力學(xué)模擬經(jīng)過一定的模擬時間達(dá)到平衡時，假設(shè)勢能被相互作用的粒子所均分，那么粒子i所具有的總能量為：ei=1

35、2mvi2+12jiij(2. 36)那么，對ei求導(dǎo)有：deidt=m vi ai+12jidijdt(2. 37)將式帶入式，整理得：Jqt=ivitei(t)+iritm vi ai+12jidijdt (2. 38)最后根據(jù)Green-Kubo線性響應(yīng)理論即可計算碳納米管在軸線方向上的導(dǎo)熱系數(shù)k。第3章 分子模擬軟件3.1 Materials Studio最強(qiáng)大的材料模擬計算平臺3.1.1 Accelrys 公司簡介Accelrys(美國)公司是世界領(lǐng)先的計算科學(xué)公司，其前身為四家世界領(lǐng)先的科學(xué)軟件公司美國Molecular Simulations Inc.(MSI)公司、Geneti

36、cs Computer Group(GCG)公司、英國Synopsys Scientific系統(tǒng)公司以及Oxford Molecular Group（OMG）公司，由這四家軟件公司于2001年6月1日合并組建的Accelrys公司是目前全球范圍內(nèi)唯一能夠提供分子模擬、材料設(shè)計以及化學(xué)信息學(xué)和生物信息學(xué)全面解決方案和相關(guān)服務(wù)的軟件供應(yīng)商，所提供的全面解決方案和科技服務(wù)滿足了當(dāng)今全球領(lǐng)先的研究和開發(fā)機(jī)構(gòu)的要求。Accelrys材料科學(xué)軟件產(chǎn)品提供了全面和完善的模擬環(huán)境，可以幫助研究者構(gòu)建、顯示和分析分子、固體、表面和界面的結(jié)構(gòu)模型，并研究、預(yù)測材料的結(jié)構(gòu)與相關(guān)性質(zhì)。Accelrys的軟件是高度模塊

37、化的集成產(chǎn)品，用戶可以自由定制、購買自己的軟件系統(tǒng)，以滿足研究工作的不同需要。Accelrys軟件用于材料科學(xué)研究的主要產(chǎn)品是Materials Studio分子模擬軟件，它可以運(yùn)行在臺式機(jī)、各類型服務(wù)器和計算集群等硬件平臺上。Materials Studio分子模擬軟件廣泛應(yīng)用在石油、化工、環(huán)境、能源、制藥、電子、食品、航空航天和汽車等工業(yè)領(lǐng)域和教育科研部門；這些領(lǐng)域中具有較大影響的跨國公司及世界著名的高校、科研院所等研究機(jī)構(gòu)幾乎都是Accelrys產(chǎn)品的用戶。3.1.2 Materials Studio的優(yōu)點（1）Materials Studio是專門為材料科學(xué)領(lǐng)域研究者開發(fā)的一款可運(yùn)行在

38、PC上的商業(yè)分子模擬軟件。支持Windows、Unix以及Linux等多種操作平臺。 （2）Materials Studio軟件采用靈活的Client-Server結(jié)構(gòu)。其核心模塊Visualizer運(yùn)行于客戶端PC，計算模塊（如Discover,Amorphous，Equilibria，DMol3，CASTEP等）運(yùn)行于服務(wù)器端，支持的系統(tǒng)包括Windows 2000、NT、SGIIRIX以及Red Hat Linux。 （3）投入成本低，易于推廣。浮動許可（Floating License）機(jī)制允許用戶將計算作業(yè)提交到網(wǎng)絡(luò)上的任何一臺服務(wù)器上，并將結(jié)果返回到客戶端進(jìn)行分析，從而最大限度地利

39、用了網(wǎng)絡(luò)資源，減少了硬件投資3.1.3 Materials Studio模塊Materials Studio采用了大家非常熟悉的Microsoft標(biāo)準(zhǔn)用戶界面，允許用戶通過各種控制面板直接對計算參數(shù)和計算結(jié)果進(jìn)行設(shè)置和分析。目前，Materials Studio軟件常用的模塊有：（1）基本環(huán)境模塊 MS.Visualizer該模塊提供了搭建分子、晶體及高分子材料結(jié)構(gòu)模型所需要的所有工具，可以操作、觀察及分析結(jié)構(gòu)模型，處理圖表、表格或文本等形式的數(shù)據(jù)，并提供軟件的基本環(huán)境和分析工具以及支持Materials Studio的其他產(chǎn)品。是Materials Studio產(chǎn)品系列的核心模塊。同時Mat

40、erials Visualizer還支持多種輸入、輸出格式，并可將動態(tài)的軌跡文件輸出成avi文件加入到Office系列產(chǎn)品中。（2）分子力學(xué)與分子動力學(xué)模塊 MS.DISCOVER MS.DISCOVER是Materials Studio的分子力學(xué)計算引擎。它使用了多種成熟的分子力學(xué)和分子動力學(xué)方法，這些方法被證明完全適應(yīng)分子設(shè)計的需要。以多個經(jīng)過仔細(xì)推導(dǎo)的力場為基礎(chǔ)，DISCOVER可以準(zhǔn)確地計算出最低能量構(gòu)象，并可給出不同系綜下體系結(jié)構(gòu)的動力學(xué)軌跡。DISCOVER還為Amorphous Cell等模塊提供了基礎(chǔ)計算方法。周期性邊界條件的引入使得它可以對固態(tài)體系進(jìn)行研究，如晶體、非晶和溶劑

41、化體系。另外，DISCOVER還提供強(qiáng)大的分析工具，可以對模擬結(jié)果進(jìn)行分析，從而得到各類結(jié)構(gòu)參數(shù)、熱力學(xué)性質(zhì)、力學(xué)性質(zhì)、動力學(xué)量以及振動強(qiáng)度。此外還有MS.COMPASS模塊， COMPASS是“Condensed-phase Optimized Molecular Potential for Atomisitic Simulation Study”的縮寫。它是一個支持對凝聚態(tài)材料進(jìn)行原子水平模擬的功能強(qiáng)大的力場。它是第一個由凝聚態(tài)性質(zhì)以及孤立分子的各種從頭算和經(jīng)驗數(shù)據(jù)等參數(shù)化并驗證的從頭算力場。使用這個力場可以在很大的溫度、壓力范圍內(nèi)精確地預(yù)測出孤立體系或凝聚態(tài)體系中各種分子的構(gòu)象、振動及熱

42、物理性質(zhì)。在COMPASS力場地最新版本中，Accelrys加入了45個以上的無機(jī)氧化物材料以及混合體系（包括有機(jī)和無機(jī)材料的界面）的一些參數(shù)，使它的應(yīng)用領(lǐng)域最終包含了大多數(shù)材料科學(xué)研究者感興趣的有機(jī)和無機(jī)材料。你可以用它來研究諸如表面、共混等非常復(fù)雜的體系。COMPASS力場是通過Discover模塊來調(diào)用的。（3）量子力學(xué)模塊 MS.DMol3和MS.CASTEPMS.DMol3是獨特的密度泛函（DFT）量子力學(xué)程序，是唯一可以模擬氣相、溶液、表面及固體等過程及性質(zhì)的商業(yè)化量子力學(xué)程序，應(yīng)用于化學(xué)、材料、化工、固體物理等許多領(lǐng)域。可用于研究均相催化、多相催化、半導(dǎo)體、分子反應(yīng)等，也可預(yù)測諸

43、如溶解度、蒸氣壓、配分函數(shù)、溶解熱、混合熱等性質(zhì)?？捎嬎隳軒ЫY(jié)構(gòu)、態(tài)密度?；趦?nèi)坐標(biāo)的算法強(qiáng)健高效，支持并行計算。在MS的高版本中加入了更方便的自旋極化設(shè)置，可用于計算磁性體系，此外還可進(jìn)行動力學(xué)計算。MS.CASTEP是進(jìn)的量子力學(xué)程序，廣泛應(yīng)用于陶瓷、半導(dǎo)體以及金屬等多種材料?？裳芯浚壕w材料的性質(zhì)（半導(dǎo)體、陶瓷、金屬、分子篩等）、表面和表面重構(gòu)的性質(zhì)、表面化學(xué)、電子結(jié)構(gòu)（能帶及態(tài)密度、聲子譜）、晶體的光學(xué)性質(zhì)、點缺陷性質(zhì)（如空位、間隙或取代摻雜）、擴(kuò)展缺陷（晶粒間界、位錯）、成分無序等?？娠@示體系的三維電荷密度及波函數(shù)、模擬STM圖像、計算電荷差分密度，高版本的MS還可以計算固體材料的紅

44、外光譜。3.2 LAMMPS強(qiáng)大的分子動力學(xué)模擬軟件LAMMPS即Large-scale Atomic/Molecular Massively Parallel Simulator，可以翻譯為大規(guī)模原子/分子并行模擬器，主要用于分子動力學(xué)相關(guān)的一些計算和模擬工作，一般來講，分子動力學(xué)所涉及到的領(lǐng)域，LAMMPS 代碼也都涉及到了。LAMMPS由美國Sandia國家實驗室開發(fā)，以GPL licence發(fā)布，即開放源代且可以免費獲取使用，這意味著使用者可以根據(jù)自己的需要自行修改源代碼。LAMMPS可以支持包括氣態(tài)，液態(tài)或者固態(tài)相形態(tài)下、各種系綜下、百萬級的原子分子體系，并提供支持多種勢函數(shù),且LA

45、MMPS 有良好的并行擴(kuò)展性。分子模擬是根據(jù)物理和化學(xué)的基本原理構(gòu)建一個模型（通常是數(shù)學(xué)模型,是對某種分子體系或反應(yīng)過程的理想化描述) ，建立一種以計算數(shù)據(jù)(由計算機(jī)來執(zhí)行)來代替實驗測量的研究方法,并獲取相關(guān)的物理和化學(xué)信息。分子模擬在材料科學(xué)方面的應(yīng)用包括模擬材料的結(jié)構(gòu)、計算材料的性質(zhì)、預(yù)測材料的行為、驗證實驗結(jié)果（重現(xiàn)實驗過程）、從微觀角度認(rèn)識材料，總之是為了更深層次理解材料的結(jié)構(gòu)，認(rèn)識材料的各種行為。分子模擬的主要優(yōu)勢在于可以降低實驗成本、具有較高的安全性、實現(xiàn)通常條件下較難或無法進(jìn)行的實驗（例如：超低溫，低于-100；超高壓，高于100Mpa）、研究極快速的反應(yīng)和變化等。 R.S.M

46、ulliken獲得諾貝爾獲獎時的感言是：“總之，我愿意強(qiáng)調(diào)我的信念：計算化學(xué)的年代已經(jīng)到來，成千上百的化學(xué)家以計算機(jī)代替實驗室，來獲得眾多的化學(xué)信息。唯一的障礙是你必須償付機(jī)時費?！?阿基米德曾說過：“給我一個支點，我就能夠翹起地球。”，但是分子模擬告訴我們:“Give me an enough powerful computer ,I can simulate the whole world”。從1980年開始，每年在Engineering Village中關(guān)于“Molecular Simulation”的文章數(shù)目由37篇遞增到最高5209（2008年）篇。與分子模擬有關(guān)的論文，美國（Uni

47、ted States）發(fā)表的篇數(shù)最多，高達(dá)16351篇，其次是日本，中國名列第三。分子模擬作為一種計算機(jī)模擬技術(shù)，主要可以進(jìn)行解釋工作和預(yù)測工作。前者為實驗奠定理論基礎(chǔ)，通過模擬解釋實驗現(xiàn)象、建立理論模型、探討過程機(jī)理等，后者為實驗過程提供可能性和可行性研究，進(jìn)行方案輔助設(shè)計、材料性能預(yù)測、過程優(yōu)化篩選等。不同的分子模擬方法可以得到不同的信息。量子力學(xué)模擬方法可以計算得到分子的大多數(shù)性質(zhì)，如結(jié)構(gòu)、構(gòu)象、偶極矩、電離能、電子親和力、電子密度、過渡態(tài)和反應(yīng)途徑等；分子力學(xué)可以計算分子體系的穩(wěn)定構(gòu)象、熱力學(xué)特性、振動光譜等；能量最小化可以探索相空間（phase space)和勢能面（pot

48、ential curve)，可以找出局部(local)與全局(global)的最小點及轉(zhuǎn)化過程的馬鞍點(saddle point)；Monte Carlo可以計算復(fù)雜分子體系的結(jié)構(gòu)變化，特別是相變化；分子動力學(xué)可以得到復(fù)雜分子的熱力學(xué)性質(zhì)、結(jié)構(gòu)、力學(xué)性質(zhì)，特別是可以觀察體系的動態(tài)演變，得到許多與時間有關(guān)的熱動力學(xué)性質(zhì)；布朗動力學(xué)可以研究蛋白質(zhì)在水溶液中的折疊過程；構(gòu)象分析可以研究復(fù)雜分子穩(wěn)態(tài)和亞穩(wěn)態(tài)結(jié)構(gòu)之間的演變等等。分子模擬所涉及的研究領(lǐng)域，涵蓋了物理、化學(xué)、化工、材料、生化等幾乎一切可以通過建立理論模型進(jìn)行研究的體系，多數(shù)能夠得到與實驗結(jié)果近似的計算結(jié)果，所以分子模擬己經(jīng)逐漸成為與實驗技術(shù)

49、并重的強(qiáng)有力的研究手段。分子模擬實際上并不僅僅局限于計算機(jī)模擬，但今天的分子模擬己和計算機(jī)模擬密不可分，正是由于計算機(jī)的高速計算技術(shù)的發(fā)展才使得分子模擬能夠像今天這樣發(fā)揮如此重要的作用。利用計算先行以了解更多的分子特性，已成為合成化學(xué)家和藥物設(shè)計學(xué)家所依賴的重要方法?；瘜W(xué)家們通過這種方法可以設(shè)計出最佳的反應(yīng)途徑，預(yù)測合成的可能性，且節(jié)省許多時間和避免材料的浪費。分子模擬除了在藥物設(shè)計方面的應(yīng)用外，已被廣泛地應(yīng)用于研究金屬材料、無機(jī)材料、高分子聚合物材料、生物材料等復(fù)雜龐大的體系。由于模擬結(jié)果隨著計算機(jī)的發(fā)展與模擬方法（算法）的改良而更加精確，分子模擬已成為如今化學(xué)各個領(lǐng)域中不可缺少的工具。1.

50、1.2 分子模擬的主要方法分子模擬的主要模擬方法有量子力學(xué)模擬和經(jīng)典力學(xué)模擬，量子力學(xué)模擬主要根據(jù)從頭算方法、半經(jīng)驗方法、DFT（Density Functional Theory,密度泛函數(shù)理論）方法，而經(jīng)典力學(xué)模擬的方法主要依據(jù)分子力學(xué)、分子動力學(xué)、蒙地卡羅模擬、布朗動力學(xué)。量子力學(xué)是利用波函數(shù)來研究微觀粒子的運(yùn)動規(guī)律的一個物理學(xué)分支學(xué)科，它以分子中電子的非定域化（delocalization）為基礎(chǔ)，一切電子的行為以其波函數(shù)（wave function）表示。根據(jù)海森伯（Heisenberg）的測不準(zhǔn)原理(uncertainly principle)，量子力學(xué)可計算區(qū)間內(nèi)電子出現(xiàn)的概率，

51、其概率正比于波函數(shù)絕對值的平方，而欲得到電子的波函數(shù)，則需解薛定諤方程式（Schrödinger equation）,即：H =E （1. 1）此式中，H為薛定諤算子（是一些數(shù)學(xué)指令）、為電子的波函數(shù)、E為能量。但是，分子中含有許多電子，要解此方程并非易事。今日，許多科學(xué)家們一方面致力于改進(jìn)量子計算的方法及增進(jìn)其精準(zhǔn)度，如1999年諾貝爾物理學(xué)獎所頒給的密度泛函數(shù)理論（density functional theory ,DFT）,即為非常精確的量子計算方法；另一方面，科學(xué)家們積極研究提升計算機(jī)的計算速度，以期計算較多電子的分子系統(tǒng)。通過量子力學(xué)，我們可以研究離子、原子、凝聚態(tài)物質(zhì)，

52、以及原子核和基本粒子的結(jié)構(gòu)、性質(zhì)，它與相對論一起構(gòu)成了現(xiàn)代物理學(xué)的理論基礎(chǔ)，也正是量子力學(xué)的出現(xiàn)，許多現(xiàn)象才得以真正被解釋，新的、無法被直覺想象出來的現(xiàn)象被預(yù)言，而且基本上所有的物理間的基本相互作用都可以在量子力學(xué)的框架內(nèi)描寫。隨著量子力學(xué)理論的逐步完善、經(jīng)驗力場的不斷開發(fā)和更快更準(zhǔn)確的算法以及計算機(jī)的普及和計算速度、計算量的不斷提升，半個世紀(jì)以來，分子模擬的理論和方法得到了快速的發(fā)展，在化學(xué)化工、材料科學(xué)、醫(yī)藥科學(xué)、生命科學(xué)等諸多領(lǐng)域發(fā)揮著越來越重要的作用，并形成了一門專門的學(xué)科分子模擬。今天，國內(nèi)外的絕大多數(shù)的優(yōu)秀科研機(jī)構(gòu)和高等院校都已或正在涉足分子模擬領(lǐng)域，每年在世界范圍內(nèi)發(fā)表的相關(guān)研究

53、論文數(shù)以萬計，并且呈不斷上升趨勢。分子模擬所擔(dān)當(dāng)?shù)慕巧灿稍缙诩兇獾慕忉屝椭饾u過渡到解釋、指導(dǎo)及預(yù)測并重型。通過上文，我們可以知道量子力學(xué)的方法適用于簡單的分子或電子數(shù)目較少的體系。但是在自然界與工業(yè)上的許多系統(tǒng)中，比如聚合物（脂肪、橡膠、安全玻璃等等）、生化分子（核酸、酶、蛋白質(zhì)、多糖類等）、合金材料等均含有大量的原子及電子。對于金屬材料、聚合物材料、濃稠溶液、固態(tài)混合物、納米材料等系統(tǒng)，我們不需要了解單一分子的性質(zhì)及分子間的相互作用，我們需要了解的是整個系統(tǒng)的各種集合的性質(zhì)、動態(tài)行為與熱力學(xué)的性質(zhì)。對于這樣的復(fù)雜體系，因其含有的電子數(shù)目較多，迄今為止，仍不可能完全仰賴量子力學(xué)計算。為解決這樣復(fù)雜龐大的體系，科學(xué)家們在二十世紀(jì)六十年代就開始著手研究各種可行的非量子力學(xué)計算方法，利用這些非量子