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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

16、ing with energy minimization）為美國加州大學的Peter Kollman等所發展，此力場通?？梢缘玫椒肿拥膸缀谓Y構、構型能、振動頻率與溶劑化自由能（Solvation free energy）等數據,所以常用于較小的蛋白質、核酸、多糖等生化分子的模擬計算。CHARMM力場（Chemistry at Harvard macromolecular mechanics）為美國哈佛大學所發展，此力場可研究許多分子系統，例如小的有機分子、溶液、聚合物、生化分子等。此力場參數除來自計算結果與實驗值的比對外，還引用了大量的量子計算結果。應用此力場能夠得到與實驗值相近的結構、作用能

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

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

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

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

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

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

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

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

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

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

27、于1931年提出。由于其解析形式簡單所以被廣泛使用，L-J作用勢的表達式為： Vrij=4rij12-(rij)6（2. 12）L-J作用勢相應的兩體作用力為：Fr=-Vr=-ddrVrijr=4(1212rij13-66rij7)r（2. 13）式中為L-J作用勢的勢能阱的參數，是L-J作用勢的互相作用勢能正好為零時的兩體距離（也稱為平衡常數）,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勢能參數，ijk表示粒子i-j鍵和粒子i-k鍵之間的夾角，鍵序bij包含多體相互作用。公式中參數的選取如表：表 2. 3 Tersoff勢函數參數eVSnm-1ancdhRnmDnm

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

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

31、的模擬時間達平衡后，分子的速度的分布近似滿足Maxwell統計速度分布，所以采取Maxwell統計分布來賦予分子的初始速度是非常理想的，Maxwell統計分布如下：pvix=(mi2kBT)1/2exp-12mivix22kBT （2. 26）該式提供了一個質量為mi的原子在溫度T時，在X方向的速度vix的高斯分布。隨著模擬時間的推進，模擬的系統將遵循能量最小化的原則達到穩定狀態。2.2.6 求解粒子的Newton運動方程在分子動力學模擬中，必須求解式的牛頓運動方程以計算速度與位置。而運動方程的迭代方法及離散方法很多，其中以Verlet算法、Gear校正預測法算法（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),所以根據式2.29，可由t及t-t的位置預測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發展了一種跳蛙方法(Leap frog method)的計算式,此方法計算速度與位置的數學表達式為：vit+12t=vit-12t+ait(2. 31)rit+t=rit+vit+12tt(2. 32)計算時假設已知vit-12t與rit，則由t時的位置rit計算質點所受的力與加速度ai(t)。再依式（）來預測時間為t+12t時的速度vit+12t，并以此類推。時間t時的速度可由式（）得出：vit=12 vit+12t+vit-12t (2. 33)利用Verlet跳蛙方法計算僅需儲存rit與vit-12t兩種

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

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

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

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

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

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

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

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

42、物理性質。在COMPASS力場地最新版本中，Accelrys加入了45個以上的無機氧化物材料以及混合體系（包括有機和無機材料的界面）的一些參數，使它的應用領域最終包含了大多數材料科學研究者感興趣的有機和無機材料。你可以用它來研究諸如表面、共混等非常復雜的體系。COMPASS力場是通過Discover模塊來調用的。（3）量子力學模塊 MS.DMol3和MS.CASTEPMS.DMol3是獨特的密度泛函（DFT）量子力學程序，是唯一可以模擬氣相、溶液、表面及固體等過程及性質的商業化量子力學程序，應用于化學、材料、化工、固體物理等許多領域?？捎糜谘芯烤啻呋?、多相催化、半導體、分子反應等，也可預測諸

43、如溶解度、蒸氣壓、配分函數、溶解熱、混合熱等性質。可計算能帶結構、態密度。基于內坐標的算法強健高效，支持并行計算。在MS的高版本中加入了更方便的自旋極化設置，可用于計算磁性體系，此外還可進行動力學計算。MS.CASTEP是進的量子力學程序，廣泛應用于陶瓷、半導體以及金屬等多種材料?？裳芯浚壕w材料的性質（半導體、陶瓷、金屬、分子篩等）、表面和表面重構的性質、表面化學、電子結構（能帶及態密度、聲子譜）、晶體的光學性質、點缺陷性質（如空位、間隙或取代摻雜）、擴展缺陷（晶粒間界、位錯）、成分無序等?？娠@示體系的三維電荷密度及波函數、模擬STM圖像、計算電荷差分密度，高版本的MS還可以計算固體材料的紅

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

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

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

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

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

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

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

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

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

53、論文數以萬計，并且呈不斷上升趨勢。分子模擬所擔當的角色也由早期純粹的解釋型逐漸過渡到解釋、指導及預測并重型。通過上文，我們可以知道量子力學的方法適用于簡單的分子或電子數目較少的體系。但是在自然界與工業上的許多系統中，比如聚合物（脂肪、橡膠、安全玻璃等等）、生化分子（核酸、酶、蛋白質、多糖類等）、合金材料等均含有大量的原子及電子。對于金屬材料、聚合物材料、濃稠溶液、固態混合物、納米材料等系統，我們不需要了解單一分子的性質及分子間的相互作用，我們需要了解的是整個系統的各種集合的性質、動態行為與熱力學的性質。對于這樣的復雜體系，因其含有的電子數目較多，迄今為止，仍不可能完全仰賴量子力學計算。為解決這樣復雜龐大的體系，科學家們在二十世紀六十年代就開始著手研究各種可行的非量子力學計算方法，利用這些非量子